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1  Introduction 
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The  broad  objective  of  the  proposed  project  is  to  investigate,  develop,  and  evaluate  computa¬ 
tionally  efficient  [1]  and  statistically  optimal  [2-4]  algorithms  for  accurate  image  reconstruc¬ 
tion  in  three-dimensional  (3D)  Ultrasonic  diffraction  tomography  (UDT)  imaging  of  the  breast 
cancer.  UDT  [5-7]  can  be  viewed  as  a  generalization  of  X-ray  tomography  where  X-rays  have 
been  replaced  with  an  acoustical  wavefield.  Because  UDT  is  non-invasive,  free  of  radiation 
hazard,  and  reproducible,  it  is  potentially  an  excellent  tool  for  imaging  of  breast  cancer  [8,9]. 
While  UDT  promises  several  potentially  important  advantages  over  conventional  ultrasonic 
imaging  and  has  found  important  uses  in  a  wide  variety  of  engineering  and  scientific  disci¬ 
plines,  its  application  to  imaging  of  breast  cancer  still  remains  largely  unexplored.  In  the  past 
two  years,  our  research  on  this  project  has  been  supported  by  a  Concept  Award  of  the  US 
Army  Medical  Research  adn  Material  Command,  and  we  believe  that  our  research  has  been 
successful  and  productive.  The  report  below  summarizes  our  research  activities  and  results  on 
the  project  to  date. 

2  Body 

Our  research  activities  on  the  project  to  date  can  be  grouped  naturally  into  4  categories.  The 
first  was  the  investigation  of  efficient  linear  algorithms  for  image  reconstruction  from  3D  data 
and  from  minimum  scan  data.  The  second  was  the  development  of  efficient  nonlinear  algo¬ 
rithms  for  2D  and  3D  image  reconstructions.  The  third  was  the  applications  of  the  developed 
algorithms  to  simulated  and  experimental  data  for  evaluation  of  their  performance.  Finally, 
as  a  by  product,  we  developed  short-scan  reconstruction  algorithms  for  reflection-mode  ultra¬ 
sound  tomography  that  can  also  be  a  potentially  important  modality  for  imaging  of  the  breast 
cancer. 

2.1  Development  of  efficient  linear  reconstruction  algorithms 

In  breast  imaging  applications  of  ultrasound,  the  first-order  Bom  or  Rytov  approximations 
are  typically  not  valid,  and  consequently,  a  linearized  Helmholtz  equation  may  not  accu¬ 
rately  describe  the  relationship  between  the  measured  scattered  wavefield  and  the  scattering 
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tissue  [10-18].  However,  many  nonlinear  reconstruction  algorithms  that  account  for  multiple¬ 
scattering  effects,  including  the  ones  we  discuss  below  [19],  involve  the  recursive  application 
of  a  linear  reconstruction  algorithm  (that  assumes  the  validity  of  the  first-order  Bom  or  Rytov 
approximation.)  It  is  therefore  very  important  to  develop  computational  efficient  and  numer¬ 
ically  robust  linear  reconstmction  algorithms  for  DT.  Below,  we  discuss  our  results  on  this 
salient  topic. 

2.1.1  Development  of  3D  efficient  reconstruction  algorithms 

In  2D  DT,  the  filtered  backpropagation  (FBPP)  algorithm  is  [5]  widely  used  for  image  recon¬ 
stmction  and  is  generally  regarded  as  being  more  accurate  than  direct  Fourier  reconstmction 
approaches.  However,  objects  such  as  the  female  breast  are  inherently  three-dimensional  and 
must  be  reconstmcted  using  fully  3D  reconstmction  algorithms  in  order  to  avoid  significant 
artifacts  and  a  loss  of  quantitative  accuracy.  We  developed  and  evaluated  novel  reconstmc¬ 
tion  algorithms  for  3D  DT,  referred  to  as  the  esstimation-combination  (E-C)  reconstmction 
algorithms,  that  effectively  solve  the  (fully)  3D  DT  reconstmction  problem  by  performing  a 
series  of  2D  Radon  transform  inversions  [20].  This  greatly  reduces  the  large  computational 
load  that  is  generally  required  by  any  other  3D  DT  reconstmction  technique  such  as  the  3D 
FBPP  algorithm.  This  is  vitally  important  for  the  development  of  computationally  tractable 
3D  nonlinear  algorithms  that  involve  the  recursive  application  of  a  3D  linear  reconstmction 
algorithm.  We  also  demonstrated  that,  in  the  presence  of  data  noise,  there  is  redundant  infor¬ 
mation  contained  in  the  3D  DT  data  function  that  can  be  exploited  by  the  3D  E-C  algorithms 
to  reduce  the  variance  of  the  reconstmcted  image  [19]. 

2.1.2  Development  of  minimum-scan  reconstruction  algorithms 

In  many  applications  of  tomographic  imaging  it  is  desirable  to  minimize  the  angular  range 
over  which  the  measurement  data  are  acquired.  This  reduces  the  time  necessary  to  collect 
the  measurement  data,  which  can  reduce  artifacts  due  to  patient  motion.  Furthermore,  in  cer¬ 
tain  situations  it  may  not  be  experimentally  possible  to  collect  data  over  a  complete  27r  range. 
We  demonstrated  that  a  minimal-scan  data  set  acquired  using  view  angles  only  in  [0,  (t>min] 
contains  all  the  information  necessary  to  reconstmct  exactly  a  2D  scattering  object  function, 
where  tt  <  ^rnin  <  37r/2  is  a  function  of  the  measurement  geometry.  Based  on  this  ob- 
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servation,  we  developed,  investigated,  and  numerically  implemented  minimal-scan  FBPP  and 
E-C  reconstruction  algorithms  for  2D  DT  that  can  exactly  reconstruct  the  scattering  object 
function  from  the  minimal  scan  data  set.  Prior  to  our  work,  all  reconstruction  algorithms 
for  DT  required  a  full  2tt  worth  of  angular  measurements  to  reconstruct  an  accurate  image. 
We  numerically  demonstrated  that  the  minimal-scan  E-C  reconstruction  algorithms  were  less 
susceptible  to  the  effects  of  data  noise  and  inconsistencies  than  were  the  minimal-scan  FBPP 
reconstruction  algorithms.  We  also  generalized  this  work  to  2D  DT  using  the  fan-beam  geom¬ 
etry  and  revealed  a  novel  relationship  between  the  maximum  scanning  angle  and  achievable 
image  resolution.  This  work  may  provide  useful  insights  into  the  development  of  minimum- 
scan  reconstruction  algorithms  for  3D  DT  that  can  be  used  for  breast  imaging  [21]. 

2.2  Development  of  efficient  nonlinear  reconstruction  algorithms 

In  certain  situations  of  the  breast  imaging,  the  first-order  Bom  or  Rytov  approximations  may 
not  be  valid.  Consequently,  a  linearized  Helmholtz  equation  may  not  accurately  describe  the 
relationship  between  the  measured  scattered  wavefield  and  the  scattering  object,  and  nonlinear 
algorithms  are  necessary  for  obtaining  accurate  images.  We  proposed  to  develop  efficient 
nonlinear  reconstmction  algorithms  for  UDT. 

2.2.1  Development  of  2D  efficient  nonlinear  reconstruction  algorithms 

Previously  we  described  our  development  and  investigation  of  E-C  reconstruction  algorithms 
for  linear  DT.  We  have  generalized  these  algorithms  to  the  case  where  a  forward  scattering 
model  includes  multiple-scattering  effects.  Two  forward  scattering  models  were  utilized  that 
captured  higher-order  terms  in  the  Bom  or  Rytov  perturbation  series,  and  are  therefore  po¬ 
tentially  useful  for  modeling  ultrasound  wave  propagation  in  breast  tissue  [22, 23].  For  each 
of  the  two  forward  scattering  models,  we  developed  families  of  nonlinear  E-C  reconstruction 
algorithms  to  solve  the  inverse  problem  [19].  The  nonlinear  E-C  reconstmction  algorithms 
operate  by  relating,  in  2D  Fourier  space  of  the  Radon  transform,  the  n-th  order  perturbation 
of  the  measured  data  function  to  the  n-th  order  perturbation  of  the  scattering  object  function. 
The  algorithms  are  recursive  in  the  sense  that  calculation  of  the  n-th  order  perturbation  of 
the  object  function  requires  knowledge  of  the  {n-l)-th  order  perturbation.  The  computational 
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efficiency  of  the  E-C  algorithms  is  therefore  very  relevant  to  this  problem.  We  also  identi¬ 
fied  consistency  conditions  for  the  nonlinear  imaging  models  employed  by  the  two  families  of 
nonlinear  E-C  algorithms.  For  both  imaging  models,  the  consistency  conditions  for  linear  DT 
were  contained  as  special  cases. 

2.2.2  Development  of  3D  efficient  nonlinear  reconstruction  algorithms 

Although  we  have  been  largely  successful  in  the  theoretical  development  of  computationally 
efficient  nonlinear  algorithms  for  2D  UDT,  the  applicability  of  such  algorithms  can  be  re¬ 
strictive  because  the  multi-scattering  effect  in  the  breast  imaging  is  generally  3D  in  nature. 
Therefore,  we  have  also  investigated  3D  nonlinear  reconstruction  algorithms.  Our  strategy  for 
the  development  of  3D  nonlinear  reconstruction  algorithms  is  similar  to  that  for  2D  nonlinear 
reconstruction  discussed  above.  Specifically,  we  proposed  to  investigate  the  two  mentioned 
forward  models  in  3D  and  to  use  the  perturbation  series  for  the  inversion.  The  inversion  of 
the  solution  at  each  perturbative  order  will  be  accomplished  through  the  use  of  our  developed 
linear  3D  E-C  algorithms  for  improving  the  computational  efficiency.  We  have  also  made 
progress  on  the  development  of  such  3D  perturbative  nonlinear  algorithms. 

2.3  Implementation  and  evaluation  of  the  proposed  algorithms 

We  have  implemented  the  proposed  linear  algorithms  and  nonlinear  algorithms  and  evaluated 
them  by  use  of  computer  simulated  data  and  real  data. 

2.3.1  Implementation  and  evaluation  of  linear  reconstruction  algorithms 

We  have  implemented  the  linear  E-C  reconstruction  algorithms  and  investigated  their  noise 
properties  by  using  a  large  number  of  computer  simulated  data  sets.  Through  our  simula¬ 
tion  studies,  we  have  demonstrated  that  it  is  possible  to  achieve  a  bias-free  reduction  of  the 
statistical  variation  in  the  reconstructed  object  function  by  utilizing  complementary  statistical 
information  inherent  in  the  scattered  data.  (The  use  of  an  explicit  smoothing  operation  gen¬ 
erally  introduces  bias  in  the  reconstructed  scattering  object  function.)  We  have  quantitatively 
demonstrated  that  the  E-C  algorithms  are  less  susceptible  to  data  noise  and  other  finite  sam¬ 
pling  effects  than  are  the  corresponding  FBPP  algorithms.  This  result  is  consistent  with  the 
observation  that  the  FBPP  algorithms  involve  more  complicated  numerical  operations  (e.g.. 
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backpropagation)  than  do  the  E-C  algorithms,  which  may  amplify  the  propagation  of  noise 
and  errors  into  the  reconstructed  image.  Using  simulated  strongly  scattering  data,  we  have 
demonstrated  that  the  E-C  algorithms  are  less  susceptible  to  modeling  errors  due  to  viola¬ 
tion  of  first-order  scattering  approximations.  These  same  results  have  been  verified  for  the 
minimum-scan  DT  problem. 

2.3.2  Implementation  and  evaluation  of  nonlinear  reconstruction  algorithms 

Using  simulated  strongly  scattering  data,  we  have  started  to  numerically  investigate  nonlinear 
reconstruction  algorithms  for  2D  DT.  As  described  in  Section  2.2.1,  our  nonlinear  reconstruc¬ 
tion  algorithms  utilize  a  forward  scattering  operator  that  assumes  the  validity  of  a  higher-order 
Bom  or  Rytov  terms.  An  accurate  numerical  implementation  of  the  forward  scattering  oper¬ 
ator  is  critical  for  obtaining  accurate  reconstructions  using  our  algorithms.  In  a  preliminary 
study,  we  have  encountered  difficulty  in  achieving  an  accurate  numerical  implementation  of 
this  operator.  The  forward  scattering  models  employed  by  our  families  of  nonlinear  algo¬ 
rithms  involve  an  integration  over  a  complex  frequency  variable,  which  is  not  computable 
in  practice.  Accordingly,  numerical  inaccuracies  were  introduced  by  truncating  the  limits  of 
integration,  which  we  observed  to  introduce  a  severe  degradation  in  the  reconstructed  image 
quality.  Compensation  for  such  image  degradations  is  a  challenging  task,  and  we  are  still 
investigating  methods  for  efficiently  and  adequately  mitigating  the  effects  of  the  integration 
tmncation  used  by  the  forward  scattering  operator. 

2.4  Development  and  evaluation  of  reconstruction  algorithms  for  reflec¬ 
tivity  tomography 

Reflectivity  tomography  has  been  applied  to  numerous  biomedical  and  non-destructive  test 
imaging  problems  [24—27].  It  has  a  strong  relationship  to  UDT  and  can  be  a  potential  useful 
technique  for  imaging  the  breast  cancer.  The  task  in  reflectivity  tomography  is  to  reconstruct 
from  the  measured  data  a  function  describing  the  reflectivity  distribution  within  the  breast. 
It  has  been  generally  considered  that  accurate  images  can  be  reconstructed  only  from  full 
scan  data  over  27r.  Recently,  we  have  investigated  and  evaluated  image  reconstruction  from 
minimum-scan  data  in  reflectivity  tomography.  Using  the  so-called  potato-peeler  perspective 
that  we  developed,  we  showed  that  accurate  images  can  be  reconstructed  from  minimum-scan 
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data  in  reflectivity  tomography.  We  also  performed  quantitative  studies  by  use  of  computer 
simulated  data,  and  the  results  in  such  studies  validated  our  theoretical  results  for  image  recon- 
stmction  in  minimum-scan  reflectivity  tomography.  One  of  our  papers  on  2D  reflectivity  to¬ 
mography  is  to  be  published  in  the  July  issue  of  IEEE  Transactions  on  Image  Processing  [28]. 
Furthermore,  we  have  generalized  our  results  to  3D  reflectivity  tomography,  and  our  paper 
on  this  topic  has  been  accepted  for  presentation  by  the  notoriously  competitive  international 
conferences  on  3D  image  reconstruction  [29]. 

3  Key  research  accomplishments 

•  We  have  developed  and  evaluated  computationally  efficient  3D  linear  reconstruction 
algorithms  that  are  more  than  100  times  faster  than  the  conventional  3D  FBPP  algorithm. 

•  We  have  investigated,  developed,  and  evaluated  algorithms  for  image  reconstruction 
from  minimum-scan  data  in  UDT  with  plane  wave  sources. 

•  We  have  investigated,  developed,  and  evaluated  algorithms  for  image  reconstruction 
from  full-scan  and  minimum-scan  data  in  UDT  with  fan-beam  wave  sources. 

•  We  have  developed  and  evaluated  computationally  efficient  3D  linear  reconstruction 
algorithms  for  UDT  with  spherical  wave  sources. 

•  We  have  developed  computationally  efficient  2D  nonlinear  reconstruction  algorithms 
for  UDT  with  plane  wave  sources. 

•  We  have  investigated  theoretically  the  development  of  computationally  efficient  nonlin¬ 
ear  reconstruction  algorithms  for  3D  UDT. 

•  We  have  developed  computer  codes  that  implement  the  proposed  linear  and  nonlinear 
reconstruction  algorithms. 

•  We  have  evaluated  the  developed  linear  and  nonlinear  reconstruction  algorithms  by  use 
of  computer  simulated  and  experimental  data. 

•  We  have  developed  and  evaluated  reconstruction  algorithms  for  2D  short-scan  reflectiv¬ 
ity  tomography. 
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•  We  have  developed  and  evaluated  reconstruction  algorithms  for  3D  short-scan  reflectiv¬ 
ity  tomography. 

4  Reportable  outcomes 

10  papers  and  3  eonference  abstracts  were  published  as  listed  in  Section  6.  Bibliography 

below. 
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data,  The  2nd  Joint  Meeting  of  the  IEEE  Engineering  in  Medicine  and  Biology  Society 
and  Biomedical  Engineering  Society,  Houston,  2002. 

3.  Y.  Zou,  M.  Anastasio,  and  X.  Pan:  Data  truncation  and  the  exterior  problem  in  reflection¬ 
mode  tomography,  IEEE  Medical  Imaging  Conference,  Norfolk,  2002. 

6  Conclusion 

Ultrasonic  diffraction  tomography  is  a  potentially  important  technique  for  imaging  of  the 
breast  cancer.  In  this  project,  we  have  investigated,  developed,  and  evaluated  computation¬ 
ally  efficient  and  statistically  optimal  algorithms  for  accurate  reconstruction  of  UDT  images 
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that  may  find  applications  in  UDT  imaging  of  breast  cancer.  In  the  past  two  years,  we  have 
made  contributions  to  UDT  research,  as  summarized  above.  Our  research  on  UDT  have  ad¬ 
dressed  numerous  scientific  and  engineering  problems  involved  in  UDT  image  reconstruction. 
These  results  are  necessary  in  making  UDT  a  viable  medical  imaging  technique  for  imaging 
breast  cancer. 
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I.  Introduction 

Because  of  its  great  potential  in  biomedical  imaging,  re¬ 
flectivity  tomography  using  ultrasound  sources  has  been 
widely  studied  since  the  late  1970s  [1-4].  Recently  there 
has  been  a  renewed  interest  in  reflectivity  tomography  due 
to  the  similarity  of  its  mathematical  structure  to  that  of  de¬ 
veloping  imaging  modalities  such  as  thermoacoustic  com¬ 
puted  tomography  (CT)  [5].  In  three-dimensional  (3D)  re¬ 
flectivity  tomography,  a  weakly  reflecting  object,  immersed 
in  an  acoustically  homogeneous  background  medium,  is  il¬ 
luminated  with  short  ultrasonic  pulses  that  are  located  at 
different  positions  on  a  spherical  surface,  enclosing  the  ob¬ 
ject,  and  the  concomitant  reflected  signals  are  measured  as 
functions  of  time  at  each  of  the  multiple  source  locations. 
The  task  in  3D  reflectivity  tomography  is  to  reconstruct 
from  such  measured  data  a  function  describing  the  reflec¬ 
tivity  distribution  within  the  scattering  object^  Under 
certain  physical  conditions,  the  3D  reflectivity  tomogra¬ 
phy  reconstruction  problem  is  tantamount  to  the  problem 
of  reconstructing  a  3D  function  from  knowledge  of  its  in¬ 
tegrals  over  sets  of  spherical  surfaces  with  varying  radii 
that  are  centered  at  the  source-receiver  locations;  i.e.,  a 
generalized  3D  Radon  transform  inversion  problem. 

In  their  seminal  paper,  Norton  and  Linzer  [2]  derived 
an  explicit  and  exact  inversion  formula  for  reconstruct¬ 
ing  the  object  function  from  data  acquired  at  all  source- 
receiver  locations  on  the  measurement  sphere.  We  refer  to 
such  an  imaging  configuration  that  utilizes  47r  steradians 
worth  of  measurements  as  a  full' scan  geometry.  In  many 
medical  imaging  applications  it  is  not  feasible  to  collect 
a  complete  data  set  using  the  full-scan  geometry.  In  this 
situation,  we  have  a  reduced-scan  measurement  geometry 
and  an  associated  limited-angle  tomographic  reconstruc¬ 
tion  problem.  When  imaging  the  female  breast,  for  ex¬ 
ample,  only  reduced-scan  measurements  acquired  over  a 
hemisphere  can  be  obtained  readily. 

In  other  tomographic  imaging  modalities  including  fan- 
beam  computed  tomography  (CT)  [6]  and  single-photon 
emission  computed  tomography  (SPECT)  [7],  it  has  been 
shown  that  one  can  accurately  reconstruct  images  from 
data  acquired  with  reduced-scan  configurations.  However, 
the  exact  inversion  formula  of  Norton  and  Linzer  [2]  cannot 

^  For  simplicity,  we  refer  to  the  object’s  reflectivity  function  as  the 
object  function. 


reconstruct,  in  general,  accurate  images  from  reduced-scan 
data.  A  reduced-scan  reconstruction  problem  for  3D  re¬ 
flectivity  tomography  has  been  investigated  in  [8].  In  that 
work,  which  was  based  on  the  paraxial  approximation,  ap¬ 
proximate  reconstruction  algorithms  were  proposed  that 
assumed  full-scan  or  reduced-scan  measurement  geome¬ 
tries.  Their  results  suggested  that  the  hemisphere  reduced- 
scan  geometry  could  yield  approximation  reconstructions 
that  were  comparable  in  quality  to  the  approximate  recon¬ 
structions  yielded  by  the  full-scan  geometry.  However,  it 
remains  unclear  whether  an  object  function  can  be  exactly 
determined  from  certain  reduced-scan  data  sets  in  3D  re¬ 
flectivity  tomography. 

In  this  work,  we  investigate  the  problem  of  reconstruct¬ 
ing  an  exact  image  from  reduced-scan  data  in  3D  reflec¬ 
tivity  tomography.  A  potato-peeler  perspective,  which  is 
conceptually  similar  to  “layer-stripping”  methods,  is  pro¬ 
posed  for  identifying  data  symmetries  in  3D  reflectivity 
tomography.  Using  the  identified  data  symmetries,  we 
heuristically  demonstrate  that,  under  certain  conditions, 
data  acquired  with  reduced-scan  configurations  in  reflec¬ 
tivity  tomography  are  sufficient  to  completely  specify  the 
object  function.  It  is  interesting  to  note  that  our  observa¬ 
tions  regarding  the  stability  of  the  reduced-scan  reflectivity 
tomography  reconstruction  problem  can  be  related  readily 
to  a  similar  analysis  of  the  limited-view  thermoacoustic  CT 
problem  [9]. 

II.  Data  Function  in  3D  Reflectivity 
Tomography 

Consider  an  acoustic  medium  that  contains  a  compactly 
supported  region  71,  residing  completely  inside  a  sphere  of 
radius  Rf  centered  at  the  origin.  The  region  R  is  charac¬ 
terized  by  an  inhomogeneous  compressibihty  /c(f)  and  den¬ 
sity  p(f).  Outside  of  R,  the  background  medium  has  no 
absorption  and  homogeneous  compressibility  kq  and  den¬ 
sity  po-  This  implies  a  constant  speed  of  sound  in  the 
background  medium  that  is  denoted  by  Cq.  The  reflectiv¬ 
ity  function  of  the  medium  is  defined  as 

/(^O  =7ife(0 -7p{f).  (1) 
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Consider  a  (measurement)  sphere,  which  has  a  radius 
Ro  and  is  centered  at  the  origin,  that  encloses  the  acous¬ 
tical  inhomogeneity,  i.e,  the  region  %.  Let  7^  €  Oq 
where  Qq  denotes  the  set  of  vectors  that  reside  on  the 
surface  of  the  measurement  sphere.  The  vector  com¬ 
ponents  of  fo  in  spherical  coordinates  are  indicated  by 
{Oq,(I>q,Ro)^  At  time  t  =  0,  an  omni-directional  acoustic 
point  source  with  time  dependence  p{l)  located  at  posi¬ 
tion  fo  emits  a  spherical  pulse  f)  into  the  region  7^, 

where  t  =  Cot  As  the  spherically-diverging  pulse  propa¬ 
gates  into  the  inhomogeneous  region  U,  echos  will  be  pro¬ 
duced  that  propagate  back  to  the  source  location.  These 
echos,  which  physically  represent  fluctuations  in  acoustic 
pressure,  are  functions  of  time  and  will  be  denoted  by  the 
function  ^(fo;f)  =  ^T(fo;^)  ”  '0tnc(fo;^>  where 
and  ^mc(fo;^  represent  the  total  and  incident  wavefields, 
respectively.  Because  V^rCfo?^  is  directly  measured  and 
'ipinc{ro;i)  is  known,  V^(fo;*)  can  be  interpreted  as  a  mea¬ 
surable  quantity. 

In  order  to  deflne  conveniently  the  relationship  between 
/(f)  and  the  measured  echo  signal  ^(fo;^),  we  introduce 
the  temporal  Fourier  transform  pairs 
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where  p  and  ^  denote  the  Fourier  transformed  quantities. 
Furthermore,  it  is  useful  to  define  the  intermediate  quan¬ 
tity 


Mro-,k)  =  ^ 
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and  let  ^i(fo;f)  denote  its  inverse  Fourier  transform.  As 
derived  in  [2],  the  mathematical  relationship  between  f{f) 
and  the  measured  echo  signal  'ip{fo,t)  is  given  by 


p(fo;?)=  f  d^ff{f)S{i-\f-fo\l  (7) 
Jn 

where  the  modified  data  function  is  g(fo;t)  =  2iipi{fo;i). 
Notice  that  g{fo;i)  is  equal  to  integrals  of  f{f)  calculated 
over  concentric  spherical  surfaces  with  radii  t  that  are  cen¬ 
tered  at  the  source-receiver  location  fo.  The  goal  of  (con¬ 
ventional  full-scan)  3D  reflectivity  tomography  is  to  utilize 
knowledge  of  p(fo;  for  fo  €  Hq  and  t  €  [0,  Ro-\-Rf]  to  de¬ 
termine  f{f}  by  inverting  the  generalized  Radon  transform 
given  in  Eq.  (7). 

Ill,  Full-scan  inversion  formula 

Let  f  =  {9,  r)  denote  a  point  inside  the  measurement 
sphere  and,  as  before,  let  fo  =  {Oo^tpo^ro)  €  Ho  reside 
on  the  measurement  surface.  The  unit  vectors  ft  cuid  no 


describe  the  directions  of  f  and  fo,  respectively.  For  a 
full-scan  geometry  where  ^(fo;  i)  is  completely  specified  for 
fo  €  Ho  and  f  €  [0,  Ro^Rf],  Norton  and  Linzer  [2]  derived 
an  exact  explicit  inversion  formula  given  by 
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In  Eq.  (8),  ;„(•),  and  P„(-)  are  the  n-th  order  spher¬ 

ical  Bessel,  Hankel  and  Legendre  polynomial  functions,  re- 
spectively. 

Due  to  the  infinite  summation,  the  numerical  implemen¬ 
tation  of  Eq.  (8)  is  computationally  demanding.  In  [2],  the 
following  approximate  inversion  formula  for  determining 
f{f)  was  also  proposed: 

/o(f)=ao/  dDo  ^(fo;2|f- fo|),  (10) 

Jqo 

where  ao  is  a  constant.  Equation  10  describes  a  simple 
backprojection  of  the  echo  data  'ip  over  the  spherical  sur¬ 
faces  from  which  the  echos  originated,  summed  over  all  of 
the  source-receiver  locations  in  Qq  (l^he  surface  of  the  mea¬ 
surement  sphere).  It  was  demonstrated  [2]  that  Eq.  (10) 
closely  approximates  the  exact  inversion  formula  given  in 
Eq.  (8).  This  result  was  explained  by  revealing  that  the 
scattering  process  itself  inherently  provides  the  correct  fil¬ 
tering  of  the  measurement  data  prior  to  application  of  the 
backprojection  operation  in  Eq.  (10). 


IV.  The  Potato-peeler  Perspective  and 
Reduced-scan  in  3D  Reflectivity 
Tomography 

The  data  function  in  Eq.  (7)  does  not  have  an  obvious 
symmetry  such  as  that  of  the  Radon  transform,  where  mea¬ 
surements  at  conjugate  views  are  mathematically  equal. 
But  if  we  assume  that  the  object  function  has  compact 
support,  then  we  can  apply  the  potato  peeler  perspective 
for  revealing  redundant  information  in  tomographic  data 
functions.  The  potato  peeler  perspective  has  been  applied 
both  heuristically  [10,11]  and  mathematically  [12,13]  to 
show  that  half-detector  [13]  and  7r-scheme  [10]  (a  gener¬ 
alization  of  the  180®  short-scan)  data  contain  enough  in¬ 
formation  to  uniquely  determine  the  object  function.  The 
same  perspective  has  also  shown  the  data  redundancy  in 
SPECT  when  the  effect  of  distant  dependent  spatial  reso¬ 
lution  is  included  [11].  The  potato  peeler  perspective  has 
already  been  applied  to  2D  reflectivity  tomography  [14]. 
We  review  the  argument  for  redund2uit  information  in  2D 
reflectivity  tomography  and  then  move  on  to  show  that 
data  covering  views  over  27r  steradians  are  sufficient  to  de¬ 
termine  the  object  function  in  3D  reflectivity  tomography. 
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Object  Function 

Fig.  1.  Schematic  of  data  function  (top)  and  two-dimensional  object 
function  (bottom)  for  the  heuristic  potato  peeler  perspective.  Points 
A,  B,  and  C  lie  on  the  outermost  ring  of  the  object  support,  and  their 
loci  in  the  data  space  are  also  shown.  The  illustrated  integration  arcs 
from  opposing  views  4>b  and  H-tt  intersect  the  support  disk  at  B. 
When  the  outermost  ring  along  with  points  A,  B  and  C  are  removed, 
D  becomes  a  point  on  the  outermost  edge. 

A.  2D  reflectivity  tomography 

In  analogy  with  the  3D  case,  the  source  position  is  spec¬ 
ified  by  a  2D  vector  in  polar  coordinates,  foicpo^Ro),  and 
Eq.  (7)  holds  as  the  data  function  except  that  the  vectors 
are  replaced  by  2D  counterparts  and  the  integration  sphere 
becomes  an  arc.  Without  loss  of  generality  we  assume  that 
the  support  of  the  object  function  lies  within  a  disk  of  ra¬ 
dius  Rf  centered  at  the  origin  of  the  measurement  circle, 
see  Fig.  1.  Let  L(^,  <f>o)  denote  the  integration  arc  from  the 
2D  reflectivity  tomography  data  function,  where  the  coor¬ 
dinate  f  -  Ho  measures  the  position  of  the  wave  fronts 
relative  to  the  center  of  the  measurement  circle.  The  key 
observation  is  that  there  exists  a  value  of  ^  such 

that  L(^Tnax,<t>B)  intersects  the  support  disk  at  one  point. 
Specifically,  consider  the  view  angle  shown  in  Fig.  1. 
The  data  point  g{fB,  -^max)  depends  on  the  object  only 
at  the  point  B;  moreover,  the  data  point  g{sB,^max)  also 


depends  only  on  B,  where  sb  —  {<Pb  +  tt,  Bo)*  Thus  the 
symmetry, 

p(r^5““4max)  ^  Oi^B'i^max)  (1-1) 

follows.  The  potato  peeler  heuristing  argument  rests  on 
this  symmetry,  and  the  fact  that  either  data  point  can  be 
used  to  determine  B. 

Conceptually,  reconstruction,  employing  the  data  redun¬ 
dancy,  proceeds  by  the  repetition  of  three  steps.  First,  each 
point  on  the  outermost  edge  of  the  disk  support  is  deter¬ 
mined  using  the  symmetry  for  the  outermost  points,  i.e. 
points  A,  B,  and  C  in  Fig.  1.  Second,  each  of  these  outer¬ 
most  points  are  forward  projected  individually  back  to  the 
data  space.  Third,  each  point’s  contribution  to  the  data  is 
subtracted  away  from  the  data  function.  When  the  sub¬ 
traction  is  performed  for  all  of  the  outermost  points,  then 
the  resulting  new  data  function  will  be  the  forward  pro¬ 
jection  corresponding  to  the  object  missing  the  outermost 
ring,  exposing  the  next  set  of  points  toward  the  interior. 
In  terms  of  Fig.  1,  once  the  contribution  of  points  A,  B, 
and  C  are  removed  from  the  data  function,  point  D  is  ex¬ 
posed  on  the  outermost  edge.  These  steps  can  be  repeated 
until  all  points  of  the  object  are  determined.  The  symme¬ 
try  used  in  the  first  step  applies  to  all  points  during  the 
peeling  procedure.  Thus,  we  have  data  redundancy,  and 
views  covering  180°  degrees  are  sufficient  to  determine  the 
object  function. 


Fig.  2.  Schematic  of  the  three-dimensional  object  function  for  the 
heuristic  potato  peeler  perspective.  The  illustrated  integration  sur¬ 
faces  from  opposing  views  fb  s^iid  5b  intersect  spherical  support  at 
one  point. 
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B.  3D  reflectivity  tomography 

The  potato  peeler  perspective  can  be  generalized  to  3D 
reflectivity  tomography  in  a  straightforward  manner.  In 
this  case  the  support  of  the  object  function  is  contained 
in  a  sphere  of  radius  Rf  shown  in  Fig.  2.  Let  the  surface 
of  integration  be  denoted  by  S(0o,  where  again  ^  = 
i-Ro  represents  the  distance  away  from  the  center  of  the 
measurement  sphere.  Clearly,  from  Fig.  2,  the  surface 
5(^o,0o,^max)  intersects  the  support  sphere  at  one  point, 
and  a  symmetry  analgous  to  Eq.  (11)  holds  for  points  on 
the  outermost  shell,  namely: 

p(^0>  =  ^(^Oj^Tnaa:)?  (^2) 

where  sb  =  (tt  —  6oj0o  +  From  this  symmetry, 

the  3D  version  of  the  heuristic  potato  peeler  perspective 
extends  naturally  from  the  2D  version.  In  the  3D  argument 
the  outermost  shell  of  the  object  function  is  determined, 
making  use  of  the  symmetry  in  Eq.  (12);  only  half  of  the 
views  over  Arr  steradians  are  needed  to  specify  all  points 
on  the  outer  shell.  Again,  the  contribution  of  points  in  the 
outer  shell  can  be  subtracted  away  from  the  data,  exposing 
the  next  shell  toward  the  interior.  From  repeating  the 
peeling,  only  27r  steradian  coverage  of  the  sphere  is  needed 
to  determine  the  complete  3D  object  function. 

V.  Numerical  Results 

We  performed  computer-simulation  studies  to  validate 
that  accurate  images  can  be  reconstructed  from  reduced- 
scan  data  in  3D  reflectivity  tomography.  A  3D  Shepp- 
Logan  phantom  was  used  to  generate  the  data  function 
For  a  given  Rq,  which  is  the  distance  between  the 
center  of  rotation  and  transducer,  g{fo;t)  is  a  function  of 
variables  where  0  <  ^  ^  0  ^  ^  27r 

specify  the  direction  of  measurements,  and  t  £  [0,  Eo  + 
Rf]  identifies  the  time  of  measurements  to  a  location.  We 
divided  the  data  space  a  60  x  120  x  128 

uniform  lattice  and  generated  full-scan  data  on  this  lattice 
for  Ro  =  30  cm. 

Since  the  data  function  in  Eq.  (7)  is  non-negative,  the  it¬ 
erative  expectation-maximization  (EM)  algorithm  can  be 
applied  to  reconstruct  the  image.  We  used  the  EM  al¬ 
gorithm  for  image  reconstruction,  because  it  is  unclear 
whether  non-iterative  algorithms  exist  for  exact  recon¬ 
struction  of  images  from  reduced-scan  data.  For  the  pur¬ 
pose  of  comparison,  we  also  used  the  EM  algorithm  for 
image  reconstruction  from  full-scan  data.  In  an  attempt 
to  increase  the  convergence  speed  of  the  EM  algorithm, 
we  used  the  ordered-subsets  EM  (OSEM)  for  image  recon¬ 
struction  from  both  full-scan  and  reduced-scan  in  our  sim¬ 
ulation  studies.  The  data  set  was  divided  into  5  subsets 
along  0  dimension,  and  the  OSEM  algorithm  was  paral¬ 
lelized  for  different  (/>  values.  We  used  200  iterations  for 
both  reconstructions  from  the  reduced-scan  and  full-scan 
data.  The  calculation  was  performed  on  the  “Chiba  City” 
cluster  computer  at  the  Mathematics  and  Computer  Sci¬ 
ence  Division  in  Argonne  National  Laboratory.  With  60 
processors,  it  took  about  80  and  40  hours  to  reconstruct 


3D  images  of  128^  voxels  from  the  full-scan  and  half-scan 
data,  respectively. 

We  show  in  Fig.  3  2D  images  at  x  =  0  cm  (left),  y  =  0  cm 
(middle),  and  z  =  -2.5  cm  (right)  in  the  reconstructed  3D 
images.  The  original  images  are  shown  in  the  top  row,  and 
the  images  in  middle  and  bottom  rows  were  obtained  from 
the  full-scan  and  half-scan  data,  respectively.  In  Fig.  4, 
we  show  the  profiles  along  the  central,  vertical  axis  of  im¬ 
ages  in  the  third  column  in  Fig.  3.  As  these  results  show, 
one  can  reconstruct  3D  images  from  half-scan  data  with 
quantitative  accuracy  comparable  to  that  of  3D  images  re¬ 
constructed  from  full-scan  data. 


Fig.  3.  2D  images  at  a;  =  0  cm  (left),  y  =  0  cm  (middle),  and 
z  =  -2.5  cm  (right)  in  the  3D  Shepp-Logan  image.  The  original 
images  are  shown  in  the  top  row.  The  images  reconstructed  from 
full-scan  and  half-scan  data  are  displayed  in  the  middle  and  bottom 
rows,  respectively. 


Fig.  4.  Original  profile  (solid  curve),  profile  in  images  obtained  from 
full-scan  (dotted  curve)  and  half-scan  (dashed  curve)  data.  Profile 
along  the  central  vertical  axes  in  images  shown  in  the  third  column 
in  Fig.  3. 
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VI.  Discussion 

We  investigated  image  reconstruction  in  3D  reflectivity 
tomography.  Motivated  by  the  potato-peeler  perspective 
that  we  developed  for  image  reconstruction  from  reduced- 
scan  data  in  2D  reflectivity  tomography,  we  proposed  a 
potato-peeler  perspective  for  identifying  symmetry  and  re¬ 
dundant  information  in  data  function  in  3D  reflectivity  to¬ 
mography.  Such  redundant  information  can  be  exploited 
for  accurate  reconstruction  3D  reflectivity  images  from 
reduced-scan  data.  We  conducted  computer  simulation 
studies,  and  results  in  these  studies  quantitatively  suggest 
that  3D  reflectivity  tomography  can  accurately  be  recon¬ 
structed  from  reduced-scan  data.  In  many  applications  of 
reflectivity  tomography,  it  is  not  feasible  to  collect  full- 
scan  data.  Therefore,  the  practical  implication  of  reflec¬ 
tivity  tomography  with  a  reduced-scan  configuration  can 
be  significant.  For  example,  when  imaging  breast  cancer, 
although  only  reduced-scan  measurements  acquired  over 
a  hemisphere  can  be  obtained,  our  results  presented  here 
suggest  that  such  data  may  contain  complete  information 
for  accurate  reconstruction  of  3D  breast  images.  We  also 
applied  the  microlocal  analysis  [15]  to  demonstrate  the  sta¬ 
bility  of  reconstructing  image  boundaries  from  reduced- 
scan  data  and  will  report  such  results  at  the  meeting.  It 
is  interesting  to  note  that  our  observations  regarding  the 
stability  of  the  reduced-scan  reflectivity  tomography  recon¬ 
struction  problem  can  be  related  readily  to  a  similar  anal¬ 
ysis  of  the  limited- view  thermoacoustic  CT  problem  [9]. 
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Full-  and  minimal-scan  reconstruction  algorithms 
for  fan-beam  diffraction  tomography 


Mark  A.  Anastasio  and  Xiaochuan  Pan 


Diffraction  tomography  (DT)  is  a  tomographic  inversion  technique  that  reconstructs  the  spatially  variant 
refractive-index  distribution  of  a  scattering  object.  In  fan-beam  DT,  the  interrogating  radiation  is  not 
a  plane  wave  but  rather  a  cylindrical  wave  front  emanating  from  a  line  source  located  a  finite  distance 
from  the  scattering  object.  We  reveal  and  examine  the  redundant  information  that  is  inherent  in  the 
fan-beam  DT  data  function.  Such  redundant  information  can  be  exploited  to  reduce  the  reconstructed 
image  variance  or,  alternatively,  to  reduce  the  angular  scanning  reqmrements  of  the  fan-beam  DT 
experiment.  We  develop  novel  filtered  backpropagation  and  estimate- combination  reconstruction  al¬ 
gorithms  for  full -scan  and  minimal-scan  fan-beam  DT.  Hie  full-scan  algorithms  utilize  measurements 
taken  over  the  angular  range  0  ^  4>  ^  2Tr,  whereas  the  minimal-scan  reconstruction  algorithms  utilize 
only  measurements  taken  over  the  angular  range  0  ^  <J)  ^  4^min»  where  tt  :s  ~  3*77/2  is  a  specified 
function  that  describes  the  fan-beam  geometry.  We  demonstrate  that  the  fiill-  and  minimal-scan  fan- 
beam  algorithms  are  mathematically  equivalent.  An  implementation  of  the  algorithms  and  numerical 
results  obtained  with  noiseless  and  noisy  simulated  data  are  presented.  ©  2001  Optical  Society  of 
America 
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1.  introduction 

Diffraction  tomography  (DT)  is  an  inversion  scheme 
that  can  be  used  for  obtaining  the  spatially  variant 
refractive-index  distribution  of  a  scattering  object. 
Applications  of  DT  can  be  found  in  various  scientific 
fields  such  as  medical  imaging, nondestructive 
evaluation  of  materials, and  geophysics.®*^  Unlike 
the  X  rays  used  in  computed  tomography  (CT),  the 
optical  or  acoustical  wave  fields  employed  in  DT  do 
not  generally  travel  along  straight  lines,  thus  pre¬ 
cluding  the  use  of  the  geometrical  optics  approxima¬ 
tion.  Therefore  a  wide  variety  of  techniques  that  are 
suitable  for  reconstruction  of  CT  images  cannot  be 
used  directly  for  reconstruction  of  diffraction  tomo¬ 
graphic  images.  CT  can  be  viewed  as  a  limiting  case 
of  DT,  in  which  the  frequency  of  the  probing  radiation 
tends  toward  infinity. 

A  vast  majority  of  the  algorithm  development  ef* 
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forts  in  DT  have  utilized  the  classic  scanning  geom¬ 
etry,*^  which  assumes  that  the  interrogating  radiation 
is  plane  wave  and  the  transmitted  scattered  wave 
field  is  measured  in  a  plane  (or  in  two  dimensions, 
along  a  line)  on  the  opposite  side  of  the  scattering 
object.  This  geometry  is  analogous  to  the  parallel- 
beam  geometry  of  x-ray  CT.  In  many  practical  sit¬ 
uations,  however,  the  interrogating  radiation  may  be 
not  plane  wave  but  rather  produced  by  a  line  source 
located  a  finite  distance  from  the  scattering  object. 
We  refer  to  this  configuration  as  the  fan-beam  geom¬ 
etry  of  DT,®  which  is  somewhat  analogous  to  the  two- 
dimensional  (2D)  fan-beam  geometry  of  CT. 

The  Bom  and  Rytov  approximations®  are  weak- 
scattering  approximations  that  effectively  linearize 
the  inhomogeneous  Helmholtz  and  Ricatti  equations, 
respectively.  The  relative  merits  of  the  Bom  and 
Rytov  approximations  in  the  context  of  DT  have  been 
widely  explored  in  the  literature.i®*^!  Under  weak- 
scattering  conditions  it  is  customary  and  useful  in  DT 
to  invoke  the  Bom  or  Rytov  approximation  that  per¬ 
mits  the  derivation  of  the  Fourier  diffraction  projec¬ 
tion  (FDP)  theorem.  The  FDP  theorem  relates  the 
one-dimensional  (ID)  Fourier  transform  of  the  mea¬ 
sured  scattered  data  to  the  2D  Fourier  transform  of 
the  scattering  object.  For  2D  DT  employing  plane- 
wave  illumination  and  the  classic  scan  configuration, 
Devaney*^  utilized  the  FDP  theorem  to  develop  the 


3334  APPLIED  OPTICS  /  Vol.  40,  No.  20/10  July  2001 


well-known  filtered  backpropagation  (FBPP)  algo¬ 
rithm,  which  can  be  viewed  as  a  generalization  of  the 
conventional  filtered  backprojection  (FBPJ)  algo¬ 
rithm  of  x-ray  CT.  Alternative  families  of  plane- 
wave  DT  reconstruction  algorithms,  referred  to  as 
estimate— combination  (E-C)  algorithms  and  general¬ 
ized  FBPP  algorithms,  have  been  developed^^  and 
investigated.!®-^'*  The  family  of  plane-wave  E-C  al¬ 
gorithms  effectively  operates  by  transforming  (in  2D 
Fourier  space)  the  DT  problem  into  a  2D  Radon 
transform  problem  that  can  be  efficiently  and  stably 
inverted  by  use  of  the  FBPJ  algorithm.  The  family 
of  plane-wave  generalized  FBPP  algorithms  recon¬ 
structs  the  image  directly  from  the  DT  data  function 
and  includes  the  FBPP  algorithm  as  a  special  mem¬ 
ber.  Both  the  generahzed  FBPP  and  E-C  algorithms 
generally  require  scattered  data  measured  from  view 
angles  in  [0, 2Tr)  to  perform  an  exact  reconstruction  of 
a  complex-valued  scattering  object.  Accordingly,  we 
refer  to  these  algorithms  as  being  fiill-scan  algo¬ 
rithms. 

Previously  we  showed*®  that,  in  plane-wave  DT 
that  employs  the  2D  classic  scanning  geometry,  a 
minimal-scan  data  set  acquired  by  use  of  view  angles 
only  in  [0,  =  3-17/2]  contains  all  the  information 

necessary  for  exact  reconstruction  of  the  scattering 
object  function.  As  the  frequency  of  the  probing  ra¬ 
diation  tends  toward  infinity,  4)mjn  ir,  which  re¬ 
flects  the  well-known  fact  that  measurements 
corresponding  to  <1)  £  [0,  -ir]  completely  specify  the  2D 
Radon  transform.  [Of  course,  compactly  supported 
objects  are  mathematically  specified  hy  a  sinogram 
p(k,  4>o).  where  <|>o  is  contained  in  any  open  set  [0,  e), 
but  ifp(5,  <t>o)  is  not  continuously  sampled  this  obser¬ 
vation  does  not  yield  stable  reconstruction  algo¬ 
rithms.]  We  subsequently  developed  minimal-scan 
FBPP*®  and  minimal-scan  E-C*®  algorithms  that 
were  capable  of  reconstructing  the  scattering  object 
function  by  use  of  the  minimal-scan  data  set.  Under 
the  conditions  of  continuous  sampling  and  in  the  ab¬ 
sence  of  noise,  we  demonstrated*®  that  the  minimal- 
scan  FBPP  and  E-C  reconstruction  algorithms  were 
mathematically  equivalent  to  the  full-scan  FBPP  and 
E-C  reconstruction  algorithms,  respectively. 

Here  we  reveal  and  examine  the  redundant  infor¬ 
mation  that  is  inherent  in  the  fan-beam  DT  data 
function.  Such  redundant  information  can  be  ex¬ 
ploited  to  reduce  the  noise  in  the  reconstructed  image 
or,  alternatively,  to  reduce  the  angular  scanning  re¬ 
quirements  of  the  fan-beam  DT  experiment.  We  de¬ 
velop  novel  E-C  and  FBPP  reconstruction  algorithms 
for  full-scan  and  minimal-scan  fan-beam  DT.  We 
demonstrate  that  the  minimal-scan  algorithms, 
which  utilize  measurements  taken  over  the  angular 
range  0  :£  <})  s  where  -ir  £  <t>mjn  -  3'**/2,  are 
mathematically  equivalent  to  their  full-scan  counter¬ 
parts  that  utilize  measurements  over  the  full  angular 
range  0  s  <J)  <  217.  An  implementation  of  the  algo¬ 
rithms  and  numerical  results  obtained  with  noiseless 
and  with  noisy  simulated  data  are  presented. 


Fig.  1.  Fan-beam  scanning  geometry  of  2D  DT.  The  interrogat¬ 
ing  cylindrical  wave  propagates  along  the  axis,  and  the  scattered 
wave  field  is  measured  along  the  line  ~  L  We  obtain  full-scan 
and  minimal-scan  data  sets  by  varying  the  measurement  angle  (j) 
between  0  and  2tt  or  between  0  and  [see  Eq.  (29)],  respectively. 
We  assume  that  S  and  D  are  much  larger  than  the  transverse 
dimension  of  the  scattering  object. 


2.  Background 

Here  we  briefly  review  the  geometry  and  approxima¬ 
tions  that  are  traditionally  employed  in  fan-beam  DT, 
as  described  in  the  pioneering  work  of  Devaney.®  A 
table  of  the  acronyms  used  in  this  manuscript  is  in¬ 
cluded  in  Appendix  A. 

A.  Fan-Beam  Diffraction  Tomography 
The  classic  scanning  configuration  of  fan-beam  DT  is 
shown  in  Fig.  1.  The  fixed  coordinate  system  {x,  y), 
the  rotated  coordinate  system  (^,  vi),  and  the  usual 
polar  coordinates  (r,  0)  are  related  by  x  =  r  cos  0,  y  - 
r  sin  0,  ^  =  X  cos  <1>  +  y  sin  <})  =  r  cos  (4)  “  0),  and  — 
-X  sin  ^  +  y  cos  <1>  =  -r  sin(<j)  -  0).  The  scattering 
object  is  illuminated  by  a  monochromatic  cylindrical- 
wave  source  located  at  the  position  iq  =  -S  on  the  -q 
axis,  emitting  a  wave  field  of  the  form 

it  TT  expO‘2'n'Vok  -  S0I) 

Uiii,  <|>)  =  l/o 

exp{j2'r7Vo[^^  -I-  (S  -I- 

(S D)T"  ’  ^ 

where  Uq  is  the^  complex  amplitude,  k  =  27rvo  is  the 
wave  number,  p  is  a  unit  vector  pointing  along  the 
positive  T]  axis,  and  D  is  the  distance  of  the  detector 
surface  from  the  center  of  rotation.  The  incident 
wave  field  4>)  could  represent  a  pressure  field  in 
acoustical  applications  or  a  scalar  electromagnetic 
field  in  optical  applications,  for  example.  From  mea¬ 
surements  of  the  scattered  wave  field  obtained  along 
the  ^  axis  oriented  at  a  measurement  angle  cj)  at  a 
distance  T[  =  D  from  the  origin,  one  seeks  to  recon¬ 
struct  the  scattering  object  function  a(r).  The  \m- 
derlying  physical  property  of  the  scattering  object 
that  is  mapped  in  DT  is  the  refractive-index  distri- 


S0I 
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bution  ra(r),  which  is  related  to  the  scattering  object 
function  by  the  equation  a(r)  =  n^(r)  -  1. 

Let  u(i,  <!))  denote  the  measured  total  wave  field 
along  the  fine  -n  =  D,  as  shown  in  Fig.  1.  The  scat¬ 
tered  data  are  given  by 

Usii,  <t>)  =  •}>)  ~  4)),  (2) 

which  can  be  treated  as  a  measurable  data  function 
because  u  (i  4))  and  u,  4))  can  be  measured.  There¬ 
fore  we  can  introduce  a  modified  data  function  M(v„, 
4)),  which  is  given  by 


M(v„,  4>) 


ITVo 


:  v'  exp[— 72tt(v'  —  Vo)0] 


X 


4>) 


‘'”'1  Ui(5.  4>) 


(3) 


where 


=  (l/2ir) 


A(?)exp(-y2Trv„5)d?, 


=  I  ^ 

^  \'S  -f-  £>’ 


(4) 

(5) 


B.  Fan-Beam  Fourier  Diffraction  Projection  Theorem 
In  plane-wave  DT,  the  FDP  theorem^'^  relates  the 
modified  data  function  to  the  2D  Fourier  transform  of 
the  scattering  object  and  can  be  viewed  as  a  general¬ 
ization  of  the  Fourier  slice  theorem  of  conventional 
x-ray  CT.  The  FDP  theorem  is  valid  under  condi¬ 
tions  of  weak  scattering  and  plane-wave  illumina¬ 
tion,  To  establish  the  FDP  theorem  for  the  fan- 
beam  DT  geometry  it  is  necessary  to  assume  the 
weak-scattering  approximation  and  the  so-called 
paraxial  approximation,®  which  is  a  well-known  ap¬ 
proximation  in  the  optics  literature.  The  paraxial 
approximation  requires  that  both  S  and  D  be  much 
larger  than  the  dimension  size  of  the  scattering  ob¬ 
ject-  This  amounts  to  requiring  that  the  largest  an¬ 
gle  subtended  by  the  object  when  the  object  is  viewed 
from  either  the  source  or  the  measurement  location 
be  much  smaller  than  a  radian. 

Under  the  Bom  and  paraxial  approximations,  Dev- 
aney  derived  the  fan-beam  FDP  theorem,®  which  is 
given  by 


Fig.  2.  The  fan-beam  FDP  theorem  states  that  <j))  is  equal 

to  the  2D  Fourier  transform  of  a(r)  along  the  semielliptic  arc  AOB 
that  has  semiaxes  equal  to  Vq  an^  ^o/x* 


As  displayed  in  Fig.  2,  Eq.  (6)  states  that  the  modified 
data  function,  M(v^,  <t)),  is  equal  to  a  semielliptical 
slice,  oriented  at  angle  (J),  through  the  2D  Fourier 
transform  of  the  object  function  a(r).  One  can  also 
derive  the  FDP  theorem  by  employing  the  Rytov  ap¬ 
proximation  instead  of  the  Bom  approximation.  In 
this  case,  Eq.  (6)  remains  unchanged  and  only  Eq.  (3) 
needs  to  be  appropriately  redefined.® 

3.  Full-Scan  Reconstruction  Algorithms  for  Fan-Beam 
Diffraction  Tomography 

First,  we  present  families  of  full-scan  FBPP  and  E-C 
reconstruction  algorithms  for  fan-beam  DT.  These 
fan-beam  algorithms  are  novel  and  contain  the  pre¬ 
viously  developed  families  of  plane-wave  FBPP  and 
E-C  algorithms  as  limiting  cases.  They  will  be  gen¬ 
eralized  to  the  minimal-scan  situation  in  Section  4 
below. 

A.  Fan-Beam  Full-Scan  Estimate-Combination 
Algorithms 

The  Radon  transform^®  of  the  scattering  object  func¬ 
tion  a(r,  6)  is  defined  as 


p(C,  4>o) 


»oc 

—  X  J  — X 


a(r,  0)8[^  -  r  cos(4)o  -  e)]dr. 


(7) 


where  4)o  is  the  projection  angle,  ^  =  r  cos(4>o  “  6).  and 
=  -r  sin(4)o  -  6).  The  2D  Fourier  transform  of p(^, 
4)o)  is  defined  as  [strictly  speaking,  P*(va)  is  the  com- 


M(v„,  <f)) 


=  j  I  o(r)exp| -y2iT  ^  ^  -  (v'  -  voIti  jdr  |v„.|  x^o. 


=  0 


|v„|  >  X  Vo¬ 


le) 
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bination  of  the  ID  Fourier  transform  with  respect  to 
Vq  and  a  ID  Fourier  series  expansion  with  respect  to 
<|)o  of  the  Radon  transform] 

X  exp[-j2irVol  -  jA(})o]d^d<|)o,  (8) 

where  is  the  spatial  frequency  with  respect  to  ^  and 
the  integer  k  is  the  angular  frequency  index  with 
respect  to  It  is  well  known  that  a{r,  6)  can 
readily  be  reconstructed  exactly  from  its  Radon 
transform  pH,  <t)o)  [or,  equivalently,  P*(yo)]  by  use  of 
a  wide  variety  of  computationally  efficient  and  nu¬ 
merically  stable  reconstruction  algorithms  such  as 
the  FBPJ  algorithm.  Therefore  the  task  of  image 
reconstruction  in  fan-beam  DT  is  tantamount  to  the 
task  of  estimating  Pa(Vo).  Furthermore,  because 
PkiVa)  for  Vq  >  0  contains  full  knowledge  of  the  Radon 
transform,  one  needs  to  estimate  P*(vo)  from  the  mea¬ 
sured  scattered  data  only  for  a  0. 


Comparison  of  Eqs.  (9)  and  (11)  yields  that,  for  |v„J 
^  XVo, 

P*(va)  =  yivJx^M.ivJ,  (13) 

provided  that 


From  Eg.  (14)  we  see  that  is  real  (that  is,  0  ^  s 
Vo  Vl  +  W)  only  for  |v,„|  <  x^o- 
In  the  absence  of  data  noise  or  inconsistencies,  one 
can  use  Eq.  (13)  to  obtain  P*(v„)  exactly  from  M^(vJ, 
which  can  readily  be  obtained  from  the  modified  data 
function.  In  the  presence  of  data  noise  or  inconsis¬ 
tencies,  one  can  use  Eq.  (13)  to  obtain  an  estimate  of 
Pkiva)-  For  any  given  0  ^  v„  ^  VqVI  +  .  we 

show  in  Appendix  B  that  four  different  roots  v„„  i  =  1, 
2,  3,  4,  satisfy  Eq.  (14).  However,  only  two  of  these 
four  roots  correspond  to  real-valued  frequencies,  which 
are  given  by 


=  -Vm2  =  V„  =  V„|  1 


1/2/ 


W)  VMl  -  (1  -  X^)(voV2vo^)  +  [1  -  X^d  -  X^)(v„Vvo^)]'^^} 


|l/2 


(15) 


From  Eqs.  (7)  and  (8)  it  can  be  shown^^  that 

Mir  fy! 

P*(v.)  =  (-»*  a(r) 

J0=O  Jr=0 

X  exp(“vA0)c/jt(2TrvQr)rdrd0,  (9) 

where  indicates  the  Ath-order  Bessel  function  of 
the  first  kind.  Because  M(v^,  (j))  is  a  periodic  func¬ 
tion  of  <J),  it  can  be  expanded  into  a  Fourier  series  with 
expansion  coefficients  given  by 


1 

MkivJ  =  —  Af(v„,  <t))exp(-jfe<|))d<|).  (10) 

2-^  Jo 


Substituting  Eq.  (6)  into  Eq.  (10),  noting  that  i  =  r 
cos(4)  -  e)  and  t|  =  -r  sin(4)  -  6),  and  carrying  out  the 
integration  over  angle  <()  (Ref.  19)  yield 


In  the  plane-wave  case  there  are  only  two  roots, 12 
which  one  can  obtain  from  Eq.  (15)  by  letting  x  ~  1- 
Therefore,  for  a  pven  0  Vo(l  +  1/x^)^^^,  one 

can  obtain  two  estimates  of  PkivJ  from  knowledge  of 
M*(v„,)  at  v„i  and  v„2,  namely. 


P*(v„) 


k 

Af*(Vml) 


7 


(16) 


P*(Va)  = 


k 

Afi(Vm2) 


M*(-vJ. 


(17) 


In  establishing  Eq.  (17)  we  used  the  property  7(“i'„') 
=  -7(v„/)“^.  It  can  readily  be  shown  that,  as  the 


M*(vJ  =  (-j)*-y( 
=  0 


I'x  |'2ir 

v„')-*  a 

Jr®0  *'8=0 


{r)exp{-jkQ)Jk(2TrryjVj„'  -  v/)rdrd0 


\Vm\  >  X^O, 


(11) 


where  vj  =  vjx^,  =  j{v'  -  Vq),  and 


'y(vm') 


vJ  + 


incident  wave  becomes  planar  (i.e.,  as  x  1)»  the 
above  results  become  the  results  in  Ref.  12  for  the 
plane-wave  case. 

In  the  absence  of  noise,  Eqs.  (16)  and  (17)  yield 
(12)  identical  (and  exact)  values  of  Pki'^a)' 

ence  of  data  noise  (or  other  inconsistencies),  the  two 
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estimates  of  PkivJ  are  distinct,  suggesting  that  it 
may  be  beneficial  to  combine  the  two  estimates  lin- 
early^^  to  form  a  final  estimate  of  P as 

P*(v„)  =  a)*(v  J[7*M*(v„)] 

+  [1  ~  c»)*(Vm)][(~l)*7  (18) 

where  is  a  combination  coefficient  that  dic¬ 

tates  how  Eqs.  (16)  and  (17)  are  combined.  This 
strategy  of  linear  combination  has  been  demon¬ 
strated  to  be  useful  in  reducing  the  noise  in  the  re¬ 
constructed  image.  12-14  Because  each  selection  of 
gives  rise  to  a  particular  final  estimate  P a(v<,), 
Eq.  (18)  can  be  interpreted  as  an  estimation  method 
for  obtaining  Pkiva)  (or,  equivalently,  the  Radon 
transform).  Because  wj(v,„)  may  be  any  complex¬ 
valued  fimction  of  v„  and  k,  Eq.  (18),  in  effect,  pro¬ 
vides  infinite  families  of  estimation  methods.  From 
the  estimate  Pki^a)  (i-O-j  Radon  transform),  one 
can  subsequently  reconstruct  the  image  a(r)  by  use  of 
the  FBPJ  algorithm.  For  simplicity,  the  use  of  Eq. 
(18)  to  estimate  Pa(Vo)  coupled  with  the  2D  FBPJ 
algorithm  to  reconstruct  a(r)  is  referred  to  as  a  fan- 
beam  full-scan  E-C  reconstruction  algorithm.  As  S 
—*  oc,  we  observe  that  x  1  nnd  that  the  fan-beam 
full-scan  E-C  reconstruction  algorithms  reduce  to  the 
plane-wave  full-scan  E-C  reconstruction  algorithms 

developed  previously.  ^2 

B.  Fan-Beam  Full-Scan  Filtered  Backpropagation 
Algorithms 

The  fan-beam  full-scan  E-C  algorithms  discussed 
above  first  estimate  Pa(v„)  (i.e.,  the  Radon  transform) 
from  the  modified  data  function  Af*(v„)  and  subse¬ 
quently  reconstruct  the  image  by  inverting  the  esti¬ 
mated  Radon  transform.  Below,  we  develop 
algorithms  that  reconstruct  images  directly  from  the 
modified  data  function.  Using  Eq.  (9),  one  can  di¬ 
rectly  express  the  object  function  in  terms  of  the  es¬ 
timate  P*(v„)  as 


a(r)  =  2-7T 


X  exp(j^6)J*(2TrVor)v„dv„. 


(19) 


Substituting  this  result  into  Eq.  (20),  and  noting  that 
v„'  =  v„/x^,  yields 

a(r)  =  ir  X  /  exp(yAe)  f  -p-;  [(1  -  X^)v' 

+  X^»'o]{(i)*(vJ-y*M*(v,„)  +  [1  -  to*(v„.)] 

X  (-1)*7-*M*(-vJ} 

X  t7*[2iT(v„'^  -  Vj,V)'''^]dv„,.  (22) 

To  reduce  Eq.  (22)  to  the  form  of  the  FBPP  algo¬ 
rithm,  we  assume  that  (Oa(v„i)  +  «»;t(  -  =  1.  Using 
the  integral  identitiesi2 

(±)*/7(v„')”*  exp(yAe)t/;fc[2-ir(v„.'^  - 


»2t7 

exp 

jk<^  ±j2lTV^-^^2TTV^J\ 

VO 

and  defining 


<t>)  =  2  exp(>A<j))<OA(vJMt(vJ 

k  =  -^ 


yield  for  Eq.  (22) 


(24) 


a(r)  =  i  f"''’  ^  [(1  -  X^)v'  +  X^o]M<“'(v,„,  4.) 
2  X  ^ 


‘'0  *'-XV0 

X  exp  j2TiVm  "2  2'irv^'ri 


d<t>dv„ 


(25) 


Equation  (25)  describes  a  family  of  fan-beam  full-scan 
FBPP  algorithms  [indexed  by  the  choice  of  <«);^(v„,)], 
which  becomes  the  family  of  FBPP  algorithms  for 
plane-wave  full-scan  DT12  when  x  In  particu¬ 

lar,  when  («),fe(v„,)  =  1/2,  Eq.  (25)  corresponds  to  the 
fan-beam  FBPP  algorithm  suggested  by  Devaney.s 
Because  the  derivation  of  the  family  of  fan-beam  full- 
scan  FBPP  algorithms  was  based  on  the  family  of 
fan-beam  full-scan  E-C  algorithms  [Eq.  (20)],  the  E-C 
and  FBPP  algorithms  are  mathematically  equiva¬ 
lent.  However,  as  will  be  demonstrated  below,  these 
FBPP  and  E-C  igorithms  respond  differently  to  data 
noise  and  other  experimental  errors. 


Substituting  Eq.  (18)  into  Eq.  (19)  yields 
*  |*vo  yJT+Th^ 

a{r)  =  2-ir  2  /  exp(y^0)  {wi(v„)7  M*(v„) 

-F  [1  -  a)*(v„,)](-l)*7'*Af*(-v,„)}JA(2Tiv„r)v„dv„. 


4.  Minimal-Scan  Reconstruction  Algorithms  for 
Fan-Beam  Diffraction  Tomography 

As  discussed  above,  the  redundant  information  con¬ 
tained  in  the  scattered  data  can  be  employed  for  re¬ 
ducing  the  image  noise.  Such  information  can  also 
be  used  for  reducing  the  angular  scanning  require¬ 
ments  in  a  DT  experiment. 


(20) 

Using  the  relationship  between  and  in  Eq.  (14), 
one  can  show  that,  for  0  :£  s  VqVI  +  1/x^ , 

v>„  =  ^^=^[(l-xV+X^o].  (21) 

V 


A.  Consistency  Conditions  and  the  Fan-Beam  Data 
Space 

According  to  the  fan-beam  FDP  theorem  in  Eq.  (6), 
the  modified  data  function  M(v^,  satisfies  the  con¬ 
sistency  condition 

Af(v^,  <t>)  =  M{-Vrn,  ^  +  TT  -  2a),  (26) 
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27r 


280 


2) 


Fig.  3.  Complete  data  space  ni"  =  U  Sft  U  U  9  contains  data 
from  the  view  angles  in  [0,  2tt\.  Subset  M  = 
obtained  from  the  view  angles  in  [0,  4)^]  is  called  the  minimal- 
complete  data  set.  The  minimal-complete  data  set  contains  all 
the  information  necessary  for  exact  reconstruction  of  the  scattering 
object  function. 


where 


sin  a  =  sgn(vm) 


(v'  -  Vo)^ 


-11/2 


(vj/xl  +  (V'  -  Vo)^J 


(27) 


Let  'M'”  =  [|v„|  ^  x^oi  0  ^  4)  s  2t7]  denote  the  complete 
(or  full-scan)  data  set.  As  shown  in  Fig.  3,  Tf  can  be 
divided  into  the  four  subspaces,  d,  and  2), 

where 


-si  =  [kml  -  Xi^o.  ,0  ^  4>  -  2a  -I-  25], 

®  =  [|v„,|  ^  x^o.  2a  +  28  <  4)  ^  -rr  -I-  2a], 
^  =  [kml  -  XVo.  -IT  +  2a  S  4.  <  4>n,in]. 

S  =  [|Vm|  ^  X^’O.  4>min  ^  4l  ^  211], 


where 


^  (1  +  i/^y72  ’  28.  (28) 

Using  Eqs.  (26)  and  (27),  one  can  verify  that  infor¬ 
mation  in  subspace  %  is  identical  to  that  in  subspace 
sd  and  that  information  in  subspace  SS  is  identical  to 
that  in  subspace  Q).  As  shown  in  Fig.  3,  because  the 
boundaries  between  the  subspaces  are  generally 
functions  of  and  <\>  and  because  each  horizontal 
line  in  V  corresponds  to  a  measurement  acquired  at 
a  particular  view  angle,  the  information  in  subspace 
^  cannot  in  practice  be  determined  independently  of 
that  in  subspace  ^  and  vice  versa.  We  therefore 
refer  to  the  union,  =  =  [|v^|  ^  ^ 

4>  —  <}>inin]j  the  minimal'Complete  data  set.  A  plot 
of  4)^in  versus  x  is  shown  in  Fig.  4.  As  x  1,  <l>min  ^ 
3tt/2,  and  M  reduces  to  the  minimal-complete  data 
set  proposed  previously  for  plane-wave  DT.^^  How¬ 
ever,  for  X  <  37r/2,  which  indicates  that  the 


Fig.  4.  Plot  of  4)inin  versus  x  reveals  that,  in  fan-beam  DT,  the 
required  angular  scanning  range  is  less  than  270°. 


angular  scanning  requirements  of  fan-beam  DT  are 
less  restrictive  than  for  plane-wave  DT. 

Figure  5  clearly  demonstrates  that  the  minimal- 
complete  data  set  contains  all  the  information  re¬ 
quired  for  an  exact  reconstruction.  According  to  the 
FDP  theorem,  the  segment  AOB  corresponds  to  a 
semielliptical  slice  through  the  2D  Fourier  transform 
of  a(r),  which  can  be  obtained  from  the  modified  data 
function  M{v^ ,  <i>) .  The  Fourier  space  coverages  pro¬ 
duced  by  the  segments  OA  and  OB  as  ^  varies  from 
0  to  are  shown  in  Figs.  5(a)  and  5(b),  respec- 


(a)  (b) 


(c) 

Fig.  5.  As  4)  varies  from  0  to  4)min»  segments  OA  and  OB 

(Fig.  2)  yield  two  incomplete  coverages  of  the  2D  Fourier  space  of 
the  object  function,  which  are  shown  in  (a)  and  (b),  respectively. 
Superimposing  the  two  incomplete  coverages  in  (a)  and  (b),  one 
obtains  a  complete  coverage,  as  shown  in  (c),  of  the  2D  Fourier 
space  of  the  object  function. 
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tively.  It  can  be  observed  that  each  of  these  two 
coverages  alone  provides  only  an  incomplete  coverage 
of  the  2D  Fourier  space  of  a(r).  However,  one  can 
superimpose  these  two  incomplete  coverages  in  Figs. 
5(a)  and  5(b)  to  obtain  complete  coverage  of  the  2D 
Fourier  space  of  a(r),  as  shown  in  Fig.  5(c). 

The  redundant  information  contained  in  subspaces 
si  and  %  of  the  minimal-complete  data  set  needs  to  be 
normalized  properly  before  or  during  the  reconstruc¬ 
tion  procedure.  Let  <{))  denote  the  minimal 

scan  data,  where  <t))  =  M{v^,  <}))  for  0  ^  ^ 

and  4))  =  0  for  <  <[.  <  2-17.  Con¬ 

sider  a  weighted  data  set  M'(v„,  <j)),  defined  as 

M'(v„,  4.)  =  u;(v„,  4>)M<'">(vm.  4)),  (29) 

where  iv{v^,  4))  can  be  a  function  of  v„  and  ^  that 
satisfies 

<|))  +  4)  +  TT  -  2a)  =  1  (30a) 

in  complete  data  space  'W* , 

U'lVm.  4>)  =  1  (3(^*^) 

in  subspace  S6,  and 

4>)  =  0  (30c) 

in  subspace  2).  One  can  choose  different  u;(v„,  4))  in 
subspaces  si  and  ^  as  long  as  these  w{v^,  4>)  satisfy 
Eq.  30(a).  In  the  numerical  examples  that  follow, 
we  used  the  explicit  form  for  w(v„,  <)))  given  by  Eq. 
(38).  We  can  now  readily  obtain  minimal-scan  re¬ 
construction  algorithms  for  fan-beam  DT. 

B.  Fan-Beam  Minimal-Scan  Estimate-Combination 
Algorithms 

Because  of  Eq.  30(c),  Eq.  (29)  can  also  be  rewritten  as 

M'(v„,  <)))  =  w{v^,  4))Af(v„,  <!>).  (31) 

Using  Eqs.  (26)  and  (30a),  one  can  verify  that 

M(v„,  4>)  =  4))  +  M'(-v„.,  4>  +  -IT  -  2a). 

(32) 

Using  Eq.  (32)  in  Eq.  (10),  one  obtains 

Mk{v„)  =  [M*'(v„)  -I-  (-1)*7"^^*'(“V,„)],  (33) 

where 


f2lT 

Mk'ivJ  =  (1/2tt)  exp{-jk^)M'{v„„  4))d4). 

Jo 

Multiplying  both  sides  of  Eq.  (33)  by  7^  and  noting 
Eq.  (16),  we  conclude  that,  for 

PkiVa)  =  [y'M.'ivJ  +  (“l)S‘"M,'(-vJ],  (34) 

where  and  are  related  by  Eq.  (14).  A  fan-beam 
minimal-scan  E-C  algorithm  is  formed  by  use  of  Eq. 
(34)  to  estimate  the  Radon  transform  from  the 
minimal-complete  data  set  acquired  at  measurement 
angles  0  <  ^  ^  <l)n,in  and  the  FBPJ  algorithm  to 
reconstruct  the  final  image.  One  can  form  different 


Fig.  6.  The  numerical  phantom  used  in  the  simulation  studies 
comprised  two  concentric  ellipses.  The  fan-beam  FDP  theorem 
was  used  to  calculate  analytically  the  simulated  scattered  field 
data  from  the  phantom. 


fan-beam  minimal-scan  E-C  algorithms  by  specif5dng 
different  choices  for  w{v„^,  <|))  in  Eq.  (31)  that  satisfy 
Eqs.  (30). 


C.  Fan-Beam  Minimal-Scan  Filtered  Backpropagation 
Algorithms 

Using  Eq.  (34)  in  Eq.  (19)  and  using  the  strategy 
outlined  in  Subsection  3.B,  we  can  also  develop  fan- 
beam  minimal-scan  FBPP  algorithms  that  can  recon¬ 
struct  the  object  function  directly  from  the  weighted 
minimal-complete  data  set  that  is  given  by 


a(r)  =  f  f  [(1  -  +  X^Vo]A^'(vm,  <l>) 

Jo  J-xvo  ^  ^ 


X  exp 


j2ttv^  ^  +  2'Trv^'Ti 


d^dv^, 


(35) 


where  ^  function  of  x  as  stated  in  Eq.  (28). 

One  can  form  different  fan-beam  minimal-scan  FBPP 
algorithms  by  specifying  different  choices  for  w(v^,  <))) 
in  Eq.  (31)  that  satisfy  Eqs.  (30).  When  x  =  1,  the 
fan-beam  minimal-scan  FBPP  algorithms  reduce  to 
the  plane-wave  minimal-scan  FBPP  algorithms  de¬ 
rived  previously.^® 


5.  Numerical  Results 

We  numerically  investigated  the  fan-beam  full-  and 
minimal-scan  reconstruction  algorithms,  using  sim¬ 
ulated  noiseless  and  noisy  data. 


A.  Data  and  Noise  Model 

We  employed  the  numerical  phantom  composed  of 
two  concentric  ellipses,  as  displayed  in  Fig.  6.  The 
values  of  the  scattering  object  function  that  corre¬ 
spond  to  the  outer  and  inner  ellipses  are  0.0005  and 
0.0001,  respectively.  We  chose  a  fan-beam  geometiy 
specified  by  x  =  0.8,  but  our  observations  below  hold 
for  arbitrary  x*  The  FDP  theorem  was  employed  to 
calculate  analytically  the  modified  data  function 
M(v^j,  4))  that  by  means  of  Eq.  (3)  determined  the 
scattered  field  data  4)).  Therefore  our  simula¬ 
tions  were  designed  to  demonstrate  the  performance 
of  the  reconstruction  algorithms  under  the  condition 
that  the  Bom  and  paraxial  approximations  are  valid. 
The  evaluation  of  the  performance  of  the  algorithms 
when  the  Bom  and  paraxial  approximations  are  not 
valid^°’^'’2^  remains  a  topic  for  future  study.  The 
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discrete  complete  data  set  comprised  128  equally 
spaced  measurement  angles  in  [0,  2it).  The  discrete 
minimal-complete  data  set  comprised  92  equally 
spaced  measurement  angles  in  [0,  4.49]  (or,  equiva¬ 
lently,  [0,  257.3°]).  In  this  way,  both  data  sets  had 
the  same  angular  sampling  increment,  A4)  =  2Tr/ 
128  «  4.49/92.  The  data  function  <j))  con¬ 

tained  129  evenly  spaced  samples  in  [-yvo,  X^oJ- 

To  simulate  the  effects  of  data  noise,  we  treated  the 
scattered  data  <}))  as  a  complex  stochastic  pro¬ 
cess  with  a  real  and  an  imaginary  component,  de¬ 
noted  u/'"’  ii,  4))  and  4)),  respectively.  (Here, 

boldface  type  for  u  and  a  denotes  a  random  variable.) 
Let  =  u/>  +.  Au/’  and  u/''  =  u +  Au 
where  uj''^  and  are  the  means  of  and  , 
respectively.  The  statistics  of  the  deviates  AUj'^ 
and  AUj^*'  are  described  by  the  circular  Gaussian 
model, 

pilLU."’,  Ak.")  - 


1 1 

X  exp 

.  +  _2 

2 

\  O'r 

/J 

(36) 


where  and  are  the  variances  of  Au^^''*  (^,  4>)  and 
Au/’(i  4)).  respectively. 

To  study  the  noise  properties  of  the  reconstructed 
images  quantitatively,  we  generated  N  =  250  noisy 
complete  and  minimal-complete  data  sets  by  using 
the  noise  model  in  Eq.  (36)  with  ~  ~  0-05.  We 

used  the  fan-beam  full-scan  and  minimal-scan  E-C 
and  FBPP  algorithms  to  reconstruct  sets  of  250  noisy 
images  from  these  noisy  data  sets.  The  matrix  size 
of  the  reconstructed  images  was  128  X  128  pixels, 
and  the  wavelength  of  the  incident  radiation  was 
equal  to  2  pixels.  The  local  image  variance  was  cal¬ 
culated  empirically  from  the  N  sets  of  reconstructed 
images  as 


var[a(r)]  = 


N-  1 


N 

i=l 


(37) 


where  ai(r)  is  the  ith  image  obtained  by  use  of  the 
reconstruction  algorithm  under  investigation. 


of Mkiv^  =  Vv„^  +  and  Ma(v„,  =  -Vv7'+^) 
from  the  sampled  values  of  and 

which  we  subsequently  used  in  Eq.  (18)  to  evaluate 
Pk{va)‘  For  each  value  of  A,  zero-padding  interpola¬ 
tion  was  employed  to  increase  the  sampling  density 
along  the  axis  of  by  a  factor  of  3  to  increase 

the  accuracy  of  the  interpolation  operation.  We  em¬ 
ployed  the  consistency  condition  Pki^a)  ” 
(-l)*P*(-vJ  to  obtain  the  64  samples  of  PkivJ  for 
spanning  the  range  VqVTTTTx  —  ^ 

implementation  of  the  FBPJ  algorithm,  an  unapo- 
dized  ramp  filter  was  used.  The  interpolation  nec¬ 
essary  for  aligning  the  backprojected  data  onto  a 
128  X  128  pixel  discrete  image  matrix  was  performed 
by  bilinear  interpolation.  When  Mf^(v^)  is  replaced 
by  Mk{v,n\  and  Eq.  (34)  is  employed  in  place  of  Eq. 
(18),  the  above  paragraph  also  describes  our  imple¬ 
mentation  of  the  fan-beam  minimal-scan  E-C  algo¬ 
rithm. 

2.  FBPP  Algorithms 

In  the  full-scan  and  minimal-scan  FBPP  algorithms, 
at  each  measurement  angle  <[),  M(v^,  <J))  or  Af'(v^,  4)), 
respectively,  was  multiplied  hy  the  depth-dependent 
filter  (|v^|/xS0[(l  ”  +  X^vo]exp[2'7rv^'n]  for  each 

of  128  discrete  values  of  Tf|.  For  each  value  of  ti,  the 
filtered  data  were  zero  padded  to  ensure  that  the 
pixel  size  of  the  reconstructed  image  matched  the 
pixel  size  of  the  images  reconstructed  by  use  of  the 
full-  and  minimal-scan  E-C  algorithms.  The  inter¬ 
polation  necessary  for  aligning  the  backpropagated 
data  onto  a  128  X  128  pixel  discrete  image  matrix 
was  performed  by  bilinear  interpolation. 

The  fan-beam  minimal-scan  E-C  and  minimal-scan 
FBPP  algorithms  utilized  the  weight  function  w{v^y 
<t))  [see  Eq.  (29)]  given  by 


4)) 

( 


sin 


2 


77  4> 

4  (7t/4)  “  CL 

77  (377/2)  4> 

4  (7t/4)  +  a 


lo 


0  <  4)  <  28  +  2a 
28  +  2a  <  4)  <  77  +  2a 

77  +  2a  <  4>  ~  4>min 

4>min  4>  ^  277 

(38) 


B.  Implementation  Details 
L  E-C  Algorithms 

From  the  uniformly  sampled  values  of  the  scattered 
field  ujiiy  4)),  Mkiv^)  can  be  determined  at  uniformly 
spaced  values  of  v^.  However,  because  of  the  non¬ 
linear  relationship  [Eq.  (14)]  between  and  v^,  the 
uniformly  spaced  values  of  at  which  is 

sampled  do  not  generally  correspond  to  the  uniformly 
spaced  values  of  at  which  one  needs  to  evaluate 
Pki^a)  utilize  the  fast  Fourier  transform  in  the 
FBPJ  algorithm).  For  each  of  65  evenly  spaced  val- 
ues  of  spanning  the  range  0  <  ^  VoVl  +  1/x^ , 

we  used  linear  interpolation  to  determine  the  values 


C.  Results 

From  the  simulated  noiseless  complete  and  minimal- 
complete  data  sets  we  reconstructed  the  phantom, 
using  the  full-  and  minimal-scan  fan-beam  recon¬ 
struction  algorithms.  Figures  7(a)  and  7(b)  show 
the  images  obtained  from  the  full-scan  and  the 
minimal-scan  E-C  reconstruction  algorithms,  respec¬ 
tively.  The  full-scan  E-C  algorithm  was  specified  by 
(jL)j^(v^)  =  1/2  in  Eq.  (18).  The  images  appear  iden¬ 
tical,  as  is  consistent  with  our  assertion  that  the  full- 
and  minimal-scan  E-C  algorithms  are  mathemati¬ 
cally  equivalent  in  the  absence  of  noise  or  other  er¬ 
rors.  Figures  7(c)  and  7(d)  show  the  images 
obtained  by  use  of  the  full-scan  and  the  minimal-scan 
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(c)  (d) 

Fig.  7.  Images  reconstructed  from  noiseless  data  by  use  of  (a) 
full-scan  E-C,  (b)  minimal-scan  E-C,  (c)  full-scan  FBPP,  and  (d) 
minimal-scan  FBPP  reconstruction  algorithms. 

FBPP  reconstruction  algorithms,  respectively. 
Again,  the  images  appear  identical,  consistent  with 
our  assertion  that  the  full-  and  minimal-scan  FBPP 
algorithms  are  mathematically  equivalent  in  the  ab¬ 
sence  of  noise  or  other  errors.  As  expected,  it  is  also 
observed  that  the  images  reconstructed  with  the  full- 
and  minimal-scan  E-C  algorithms  [Figs.  7(a)  and 
7(b)]  are  identical  to  the  images  reconstructed  with 
the  full-  and  minimal-scan  FBPP  algorithms  [Figs. 
7(c)  and  7(d)]. 

Using  one  of  the  simulated  noisy  complete  and 
minimal-complete  data  sets,  we  again  reconstructed 
the  phantom,  using  the  full-  and  minimal-scan  E-C 
and  FBPP  reconstruction  algorithms.  Figures  8(a) 
and  8(b)  show  the  images  obtained  by  use  of  the 
full-scan  and  the  minimal-scan  E-C  reconstruction 
algorithms,  respectively.  The  images  no  longer  ap¬ 
pear  identical,  and  the  image  reconstructed  with  the 
full-scan  E-C  algorithm  appears  less  noisy  than  the 
image  reconstructed  with  the  minimal-scan  E-C  al¬ 
gorithm.  Figures  8(c)  and  8(d)  show  the  images  ob¬ 
tained  by  use  of  the  full-scan  and  the  minimal-scan 
FBPP  reconstruction  algorithms,  respectively.  Simi¬ 
larly,  the  image  reconstructed  with  the  full-scan 
FBPP  algorithm  appears  less  noisy  than  the  image 
reconstructed  with  the  minimal-scan  FBPP  algo¬ 
rithm. 

The  observation  that  the  full-scan  algorithms  gen¬ 
erate  cleaner-looking  images  than  do  the  minimal- 
scan  algorithms  is  not  surprising  and  can  be 
qualitatively  understood  by  examination  of  ways  in 
which  the  redundant  information  inherent  in  the  DT 
data  function  is  utilized.  The  full-scan  FBPP  and 
E-C  algorithms  [with  ^^(v^)  7^  0,  1]  effectively  use 
the  redundant  information  to  reconstruct  two  sepa- 


(c)  (d) 

Fig.  8.  Images  reconstructed  from  noisy  data  by  use  of  (a)  full- 
scan  E-C,  (b)  minimal-scan  E-C,  (c)  full-scan  FBPP,  and  (d) 
minimal-scan  FBPP  reconstruction  algorithms.  The  noisy  data 
were  generated  with  cr,.  =  o,  =  0.05  in  Eq.  (36). 

rate  images  that  are  averaged  to  form  the  final  image. 
It  has  been  quantitatively  demonstrated  that  this 
effective  averaging  operation  can  result  in  an  unbi¬ 
ased  reduction  of  the  reconstructed  image  vari¬ 
ance.  12.14  The  minimal-scan  algorithms,  however, 
utilize  part  of  the  redundant  information  that  is  in¬ 
herent  in  the  data  function  to  reduce  the  angular 
range  over  which  measurements  are  required  for  the 
reconstruction.  The  redundant  information  not 
used  for  this  purpose  can  be  used  to  reduce  the  image 
variance  of  the  reconstructed  image.  Specifically, 
the  complementary  information  contained  in  sub¬ 
spaces  si  and  %  of  Fig.  3  is  weighted,  as  described  by 
Eq.  (29),  and  is  subsequently  combined  during  the 
reconstruction  procedure.  However,  unlike  the  full- 
scan  algorithms,  the  minimal-scan  algorithms  cannot 
further  reduce  the  reconstructed  image  variance  by 
exploiting  the  fact  that  subspaces  ^  and  Q)  contain 
redundant  information. 

Although  the  full-  and  minimal-scan  E-C  algo¬ 
rithms  are  mathematically  equivalent  to  the  full-  and 
minimal-scan  FBPP  algorithms,  we  observed  from 
Fig.  8  that  the  E-C  and  FBPP  algorithms  respond 
differently  to  noise  that  is  present  in  a  discrete  data 
set.  To  confirm  this  observation  quantitatively,  we 
calculated  the  local  image  variances  of  images  recon¬ 
structed,  using  the  different  methods.  Figure  9(a)  is 
a  plot  of  the  local  variance  obtained  from  the 
minimal-scan  FBPP  reconstructed  images  divided  by 
the  local  variance  obtained  from  the  minimal-scan 
E-C  reconstructed  images.  Clearly,  the  ratio  of  the 
variances  is  ever3where  greater  than  1,  quantita¬ 
tively  demonstrating  that  the  minimal-scan  E-C  re¬ 
construction  algorithms  are  less  susceptible  to  the 
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(b) 

Fig.  9.  (a)  Plot  of  the  local  variance  obtained  from  the  full-scan 
FBPP  reconstructed  images  divided  by  the  local  variance  obtained 
from  full-scan  E-C  reconstructed  images,  (b)  Plot  of  the  local 
variance  obtained  from  the  minimal-scan  FBPP  reconstructed  im¬ 
ages  divided  by  the  local  variance  obtained  from  minimal-scan  E-C 
reconstructed  images.  In  both  the  full-  and  the  minimal-scan 
cases,  the  E-C  algorithms  did  a  better  job  of  suppressing  data  noise 
than  did  the  FBPP  algorithms. 


effects  of  data  noise  than  are  the  minimal-scan  FBPP 
reconstruction  algorithms.  Figure  9(b)  is  a  plot  of 
the  local  variance  obtained  from  the  full-scan  FBPP 
reconstructed  images  divided  by  the  local  variance 
obtained  from  the  full-scan  E-C  reconstructed  im¬ 
ages.  The  ratio  of  the  variances  is  everywhere 
greater  than  1,  quantitatively  demonstrating  that 
the  full-scan  E-C  reconstruction  algorithms  are  less 
susceptible  to  the  effects  of  data  noise  than  are  the 
full-scan  FBPP  reconstruction  algorithms.  These 
results  are  consistent  with  those  obtained  previously 
for  full-  and  minimal-scan  DT  utilizing  plane-wave 

illumination. 

6.  Summary 

In  this  study,  we  revealed  and  examined  the  redun¬ 
dant  information  that  is  inherent  in  the  fan-beam  DT 


data  function.  Such  information  can  be  exploited  to 
reduce  the  reconstructed  image  variance  or  alterna¬ 
tively  to  reduce  the  angular  scanning  requirements  of 
the  fan-beam  DT  experiment.  We  developed  novel 
full-scan  and  minimal-scan  E-C  and  FBPP  recon¬ 
struction  algorithms  for  fan-beam  DT.  The  family  of 
fan-beam  full-scan  E-C  algorithms  operates  by  trans¬ 
forming  (in  2D  Fourier  space)  the  fan-beam  DT  prob¬ 
lem  into  a  2D  parallel-beam  x-ray  CT  problem,  which 
can  be  efficiently  and  stably  inverted  by  use  of  the 
FBPJ  algorithm.  The  family  of  fan-beam  full-scan 
FBPP  algorithms  operates  directly  on  the  modified 
data  function  to  reconstruct  the  image  and  contains 
the  fan-beam  FBPP  algorithm  suggested  by  Dev- 
aney®  as  a  special  member.  Different  members  of 
the  families  of  full-scan  E-C  and  FBPP  algorithms 
are  specified  by  different  choices  of  the  combination 
coefficient  which  controls  ways  in  which  the 

redundant  information  in  the  data  ffinction  is  com¬ 
bined.  Reconstruction  algorithms  that  correspond 
to  different  choices  of  will  in  general  respond 

differently  to  the  effect  of  noise  and  discrete  sam¬ 
pling. 

The  fan-beam  minimal-scan  E-C  and  FBPP  algo¬ 
rithms  were  developed  from  the  concept  of  the 
minimal-complete  data  set.  The  minimal-complete 
data  set,  which  is  acquired  by  use  of  view  angles  only 
in  [0,  where  tt  <  ^^nin  ^  37r/2,  contains  all  the 
information  necessary  for  exactly  reconstructing  the 
scattering  object  function.  The  fan-beam  minimal- 
scan  E-C  and  FBPP  algorithms  utilize  a  weighting 
function  w{v^,  <t>)  to  normalize  appropriately  the  par¬ 
tially  redundant  information  inherent  in  the 
minimal-complete  data  set.  Accordingly,  one  can 
form  different  fan-beam  minimal-scan  E-C  and  FBPP 
algorithms  by  specifying  different  choices  for  this 
weighting  function.  Reconstruction  algorithms  that 
correspond  to  different  choices  of  w{v^,  <{))  will  in 
general  respond  differently  to  the  effect  of  noise  and 
discrete  sampling.^®  It  can  be  readily  verified  that, 
under  the  conditions  of  continuous  sampling  and  in 
the  absence  of  noise,  the  minimal-scan  E-C  and  FBPP 
algorithms  are  exact  and  mathematically  equivalent 
to  their  full-scan  counterparts  that  utilize  measure¬ 
ments  over  the  angular  range  0^4)^  2tt. 

An  implementation  of  the  fan-beam  full-scan  and 
minimal-scan  algorithms  has  been  presented,  along 
with  numerical  results  obtained  with  noiseless  and 
with  noisy  simulated  data.  It  was  observed  that  the 
full-scan  algorithms  did  a  better  job  of  suppressing 
data  noise  than  did  their  minimal-scan  counterparts. 
We  quantitatively  demonstrated  that  the  full-  and 
minimal-scan  E-C  algorithms  are  less  susceptible  to 
data  noise  and  to  finite  sampling  effects  than  are  the 
full-  and  minimal-scan  FBPP  algorithms,  respec¬ 
tively.  This  result  is  consistent  with  the  observation 
that  the  FBPP-based  algorithms  involve  more- 
complicated  numerical  operations  than  do  the  E-C- 
based  algorithms,  which  may  amplify  the 
propagation  of  noise  and  errors  into  the  reconstructed 
image. 

We  have  assumed  a  2D  imaging  model  in  this 
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study.  Therefore  the  developed  reconstruction  algo¬ 
rithms  may  be  useful  for  applications  in  which  out- 


It  can  readily  be  verified  that  the  four  roots  of  Eq. 
(A2)  are  given  by 


V;„l-  V^2”VJ1  J  -  (1  ^  x')(VaV2Vo")  +  [1  '  X'd  - 


1/2 


Vm3 


[V2{1  -  (1  -  x')(VaV2vo')  -  [1  -  x'(l  -  x')(v.Vvo')]'^'}, 


1/2 


(A3) 

(A4) 


of-plane  scattering  is  not  significant.  The  full-scan 
E-C  and  FBPP  reconstruction  algorithms  can  be  gen¬ 
eralized  readily  to  address  the  three-dimensional  DT 
problem  by  use  of  spherical-wave  sources  and  planar 
measurement  surfaces.  It  remains  unclear  whether 
numerically  stable  versions  of  the  minimal-scan  E-C 
and  FBPP  reconstruction  algorithms  can  be  devel¬ 
oped  for  three-dimensional  imaging  geometries. 

Here  we  have  developed  linear  reconstruction  al¬ 
gorithms  for  fan-beam  DT.  It  was  not  our  intent  to 
address  the  limitations  of  the  Bom  or  Rytov^®'^^ 
weak-scattering  approximation.  The  developed  full- 
and  minimal-scan  algorithms  will,  however,  provide 
a  natural  framework  for  the  incorporation  of  higher- 
order  scattering  perturbation  approximations22-24 
into  the  algorithms.  It  remains  to  be  determined 
whether  minimal-scan  reconstruction  algorithms  can 
be  developed  without  use  of  the  paraxial  approxima¬ 
tion. We  intend  to  report  on  the  theoretical  devel¬ 
opment  and  numerical  analysis  of  these  problems  in 
a  forthcoming  publication. 


Appendix  A:  Acronyms  Used 


CT 

DT 

E-C  reconstruction 
algorithm 

FBPJ  reconstruction 
algorithm 

FBPP  reconstruction 
algorithm 

FDP  theorem 


Computed  tomography 
Diffraction  tomography 
Estimate-combination  recon¬ 
struction  algorithm 
Filtered  backprojection  re¬ 
construction  algorithm 
Filtered  backpropagation  re¬ 
construction  algorithm 
Fourier  diffraction  projection 
theorem 


Appendix  B:  Relationships  between  and 

From  Eq.  (14)  we  know  that 


For  a  given  that  satisfies  0  <  <  Vq  Vl+Vx^ » we 
would  like  to  find  the  values  of  that  satisfy  Eq. 
(Al) .  This  is  equivalent  to  solving  for  the  roots  of  the 
fourth-order  equation: 

(1  -  +  [W  -  2(1  -  xW](^)' 

+  v<.W-W)  =  0.  (A2) 


Because  when  0  ^  ^  VqVI  +  1/x^ 


1-(1- 


1  -  x"(i  -  X^)  ^ 

Vo 


1/2 


(A5) 


we  observe  that  the  roots  v„3  and  v„4  are  complex 
valued  and  therefore  are  not  physically  meaningful. 
As  expected,  when  x  1>  EQ-  (A3)  reduces  to  the 
known  result^^  for  plane-wave  DT  given  by 

v„i=-v„2  =  v„|l-^j  .  (A6) 
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An  Improved  Reconstruction  Algorithm  for  3-D  Diffraction 
Tomography  Using  Spherical- Wave  Sources 

Mark  A.  Anastasio*  and  Xiaochuan  Pan 


Diffraction  tomography  (DT)  is  an  inversion  technique  that 
reconstructs  the  refractive  index  distribution  of  a  weakly  scattering  object. 
In  this  paper,  a  novel  reconstruction  algorithm  for  three-dimensional 
diffraction  tomography  employing  spherical-wave  sources  is  mathe¬ 
matically  developed  and  numerically  implemented.  Our  algorithm  is 
numerically  robust  and  is  much  more  computationally  efficient  than 
the  conventional  filtered  backpropagation  algorithm.  Our  previously 
developed  algorithm  for  DT  using  plane-wave  sources  is  contained  as 
a  special  case. 

Index  Diffraction  tomography,  wavefield  inversion,  3-D  imaging. 


I.  Introduction 

In  diffraction  tomography  (DT),  a  weakly  scattering  object  is  interro¬ 
gated  using  a  diffracting  wavefield,  and  the  scattered  wavefield  around 
the  object  is  measured  and  used  to  reconstruct  the  (low-pass  filtered) 
refractive  index  distribution  of  the  scattering  object.  The  principles  of 
DT  have  been  extensively  utilized  for  developing  optical  [I],  [2]  and 
acoustic  [3]  tomographic  imaging  systems  for  biomedical  applications. 

It  is  known  that  the  filtered  backpropagation  (FBPP)  and  direct 
Fourier  (DF)  reconstruction  algorithms  for  three-dimensional  (3-D) 
DT  possess  certain  limitations  [4].  The  depth-dependent  filtering 
(backpropagation)  in  the  3-D  FBPP  algorithm  requires  a  large  number 
of  two-dimensional  (2-D)  fast  Fourier  transforms  (FFTs)  to  be  per¬ 
formed  for  processing  the  measured  data  at  each  measurement  view, 
which  renders  the  3-D  FBPP  algorithm  computationally  demanding. 
Furthermore,  we  have  shown  that  [in  two-dimensional  (2-D)  DT] 
the  FBPP  algorithm  may  amplify  data  noise  more  than  alternative 
algorithms  do  [5].  The  3-D  DF  algorithms  require  the  use  of  a  3-D 
interpolation  method  to  obtain  samples  on  a  3-D  Cartesian  grid  in 
the  Fourier  space  of  the  scattering  object,  upon  which  a  3-D  inverse 
FFT  can  be  employed  to  reconstruct  the  scattering  object  function. 
Because  the  sample  density  in  the  3-D  Fourier  space  obtained  from 
the  measured  data  is  nonuniform,  sophisticated  and  computationally 
demanding  interpolation  strategies  are  generally  required  to  avoid 
producing  significant  interpolation  errors  that  would  degrade  the 
accuracy  of  the  reconstructed  image. 

For  3-D  DT  employing  plane-wave  sources,  we  have  recently 
developed  a  new  class  of  reconstruction  algorithms  that  circumvent 
the  shortcomings  of  the  3-D  FBPP  and  DF  algorithms  [4].  These 
algorithms,  referred  to  as  plane-wave  estimate-combination  (E-C)  re¬ 
construction  algorithms,  effectively  reduce  the  3-D  DT  reconstruction 
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problem  to  a  series  of  2-D  X-ray  reconstruction  problems  and,  thus, 
greatly  reduce  the  large  computational  load  required  by  conventional 
3-D  DT  reconstruction  algorithms.  Additionally,  these  algorithms  do 
not  require  an  explicit  3-D  interpolation  in  the  Fourier  space  of  the 
scattering  object. 

In  many  imaging  applications  [1],  [6],  it  may  be  useful  to  utilize  a 
diverging  spherical-wave  rather  than  a  plane-wave  to  interrogate  the 
scattering  object.  Because  of  the  distinct  advantages  of  the  E-C  recon¬ 
struction  algorithms  for  plane-wave  DT  [4],  it  is  important  to  gener¬ 
alize  them  to  DT  employing  spherical-wave  sources.  In  this  paper,  we 
generalize  our  previously  developed  (plane-wave)  E-C  reconstruction 
algorithms  to  DT  employing  spherical-wave  sources  and  numerically 
demonstrate  the  developed  algorithm  using  simulated  data. 


II.  Background 


We  will  utilize  the  model  of  spherical-wave  DT  described  by  De- 
vaney  in  [7],  the  scanning  geometry  of  which  is  shown  in  Fig.  1 .  The 
case  of  3-D  DT  utilizing  a  plane-wave  source  can  be  viewed  as  a  special 
case  of  spherical-wave  DT  and  will  be  discussed  below.  The  scattering 
object  is  illuminated  by  monochromatic  spherical-wave  source  located 
at  the  position  ??  =  -  S  on  the  r/  axis,  emitting  a  wavefield  of  the  form 


-2,  <!>)  —Ao 


^j2wuo\r-S^\ 

\r-  5/j| 


(1) 


where  Ao  is  the  complex  amplitude,  k  =  27r^'o  is  the  wavenumber, 
is  a  unit  vector  pointing  along  the  positive  r)  axis,  and  D  is  the  distance 
of  the  detector  surface  (i.e.,  i~z  plane)  to  the  axis  of  rotation  (i.e., 
axis).  A  complete  data  set  is  acquired  by  varying  the  view  angle  6  be¬ 
tween  0  and  27r,  where  <!>  denotes  the  angle  measured  from  the  positive 
X  axis.  (The  rotated  coordinates  (^ ,  ?/)  are  related  to  the  unrotated  co¬ 
ordinates  (x,  1/)  by  C  =  ^  +  2/  sin  6  and  ?/  =  2/  cos  0  -  ar  sin  <6.) 

From  the  measured  scattered  wavefield  data,  one  seeks  to  reconstruct 
the  scattering  object  function  a(f ).  The  underlying  physical  property 
of  the  scattering  object  that  is  mapped  in  DT  is  the  refractive  index  dis¬ 
tribution  n  (f* ),  which  is  related  to  the  scattering  object  function  by  the 
equation  fl.{f )  =  n^(r )  —  1. 

Let  0)  denote  the  measured  total  wavefield  in  the  i-z  plane 
positioned  at  =  D,  as  shown  in  Fig.  1 .  The  scattered  data  is  given  by 


Us{i,z,<l>) —u{i,z,6)  -  uq{^,z,<I>)  (2) 


which  can  be  treated  as  a  measurable  data  function  because  «(^,  i:,  0) 
can  be  measured  and  because  2,0)  is  assumed  known.  We  can 
introduce  a  modified  data  function 


tt’S  2,  0) 
lio(^,2,0) 


} 


(3) 


where:r..„,  Mi,  ^)}  =  (l/2>r)  d?  dz 


and 


X  = 


/ 


5 

S  +  D 


(4) 


(5) 
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Fig.  1.  The  classical  scanning  geometry  of  spherical -wave  3-D  DT.  A 
spherical-wave  is  incident  at  angle  <i>,  and  the  scattered  field  is  measured  in 
the  plane  positioned  at  r/  =  JD.  The  measurement  angle  is  varied  from 

0to27r. 


Equation  (6),  coupled  with  a  3-D  interpolation  method,  can  be  used 
to  implement  a  DF  reconstruction  algorithm.  The  3-D  FBPP  algorithm 
for  plane-wave  DT  [8]  can  readily  be  extended  to  3-D  spherical-wave 
DT  [6]  and  is  given  by  [9],  [6] 

X  M(^„,  r..,  dp„.  dv..  (7) 

where 


and  _ _ 

V, = -j  . 

III.  E-C  Reconstruction  Algorithm 

In  deriving  the  E-C  reconstructions  algorithms  for  spherical-wave 
DT,  we  will  modify  the  procedure  employed  for  the  derivation 
of  the  plane- wave  DT  algorithms  described  in  [4].  Let  a(r,  0.  c) 
denote  a  3-D  scattering  object  function  where  r  =  -f 

B  =  tan“^(2//a:),  and  ^  denote  cylindrical  coordinates.  The  X-ray 
transform,  p(^,  -s,  do),  of  a(r,  0,  z)  is  defined  as 


pU,  ^7  do) 


-/■ 

Jt}~- 


a(r,  S,  z)d7} 


(8) 


where  do  is  the  projection  angle,  i  =  rcos(0o  -*  ^)  and 
r]  =  ~rsiii(0o  “  ^)-  Equation  (8)  states  that  p(^,  0o)  is 

the  geometrical  projection  of  a(r,  z)  onto  the  (~z  plane  oriented 
at  angle  do  about  the  2  axis.  Consequently,  p(^,  2,  do )  can  be 
interpreted  as  a  stack  of  2-D  Radon  transforms  of  .c)  along  the 

axis.  The  combination  of  a  2-D  Fourier  transform  with  respect  to  ^ 
and  z  and  a  one-dimensional  Fourier  series  expansion  with  respect  to 
do  of  p(C7  <^7  do)  is  given  by 


Fig.  2.  A  plot  of  2/m  as  a  function  of  and  i^a  as  described  by  (17).  The  plot 
was  generated  using  \  =  0.7  and  uo  =  t*. 


r /r  p(c.^,^o) 

J^,z=~-cx> 


Using  paraxiaP  and  weak  scattering  approximations,  Devaney  de¬ 
rived  [7]  the  spherical-wave  FDP  theorem  that  is  given  by 

d) 

./  ~  3C  J  —  X 

=  0  ifi^m  +  (^) 


where  the  integer  k  is  the  angular  frequency  index  with  respect  to  do, 
andand  are  the  continuous  spatial  frequencies  conjugate  to  C  and 

z,  respectively.  As  a  matter  of  convenience,  we  refer  to  Pfc(i^a,  I'z)  as 
the  “5-D  Fourier  transfornF  of  p(^,  z,  do).  Substituting  (8)  into  (9) 
and  carrying  out  the  integral  over  do  yields 

P.K,  ..)  =  {-if  r  c— dz  r  r  a(r,  z) 

J z~  —  oc  JQ—Q  j r—Q 

xe•'•'^^]Jt{27^l/„r)rdrd^.  (10) 


Equation  (6)  states  that  the  modified  data  function,  M(j/rrM  ^2,  <P)y 
is  equal  to  a  semi-ellipsoidal  slice,  oriented  at  angle  d,  through  the 
3-D  Fourier  transform  of  the  object  function  a(f ).  As  the  measure¬ 
ment  angle  d  is  varied  between  0  and  27r,  a  certain  low-pass  coverage 
of  the  3-D  Fourier  space  is  specified  (for  details,  please  refer  to  refer¬ 
ence  [7]).  In  this  paper,  we  assume  the  validity  of  the  first-order  Bom 
approximation;  the  spherical -wave  FDP  theorem  can  also  be  derived 
based  upon  the  Rytov  approximation  by  appropriately  redefining  (3) 
[7]. 

*The  paraxial  approximation  requires  that  both  5  and  D  are  much  larger  than 
the  size  of  the  scattering  object. 


Again,  for  convenience,  we  refer  to  1^2)  =  (1/^^) 

M(i^m,  Vz,  dd  as  the  “i-D  Fourier  transform**  of  the 

modified  data  function.  Using  (6)  and  noting  that  $  =  r  cos(d  - 
and  7}  =  -r  sin(d  —  one  can  obtain 


I/-) 

/oc  /*2t  roc 

/  /  a{r,e,z) 

=  -oc  J9-0  ./r=0 

X  I  ^  ^.2.rt.^rsm(<t.-<»)->(A:*+2»..;„rcos(*-(>)) 


(11) 
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where  •  Carrying  out  the  integral  [10]  in  the  curly 

brackets  in  ( 1 1 ),  one  can  re-express  M/c  (vm  <,  r'z )  as 


(12) 


where  v'i,  + 1/?  <  x^i'o  ■  ,  ■>  2  2 

Comparison  of(lO)  and  (12)  yields  that  for  i/„ +1^..  <  x  vo 

Pki-'a.  =  bVin,  V..)  (13) 


where 


The  explicit  forms  of  this  general  relationship  in  (13)  are  determined 
completely  by  the  forms  of  the  relationship  between  i^a  and  Um  implied 
in  (14).  Equation  (14)  indicates  that  the  condition  i/l,  +  <  X't'o 

is  equivalent  to  the  condition  i^a  +  ^  i^o  (1  +  order  to 

determine  as  a  function  of  i^z  and  we  must  solve  for  the  roots 
of  the  fourth-order  equation 


+C2i''J  +  C3  =  Q  (16) 


where  the  coefficients  C,  are  given  by 

C,  =(l-x^)^ 

C2  =2(1  -  x^)  (21/i;  -  u'i  -  +41-? x'" 


and  also  from  I'r)  at  i'm2  =  —I'm  and  as 

Pfc(l'a,  ’'•)  =[7(''".2,  I'-) 

=  (-l)*[7(«'m,  »'.)•  (21) 

The  fact  that  there  are  two  distinct  ways  to  obtain  (i^a .  t'z)  from  the 
measured  data  can  be  explained  by  the  existence  of  a  double  coverage 
of  the  scattering  object  Fourier  space  [7],  At  a  given  measurement  angle 
d,  the  spherical-wave  FDP  theorem  relates  the  2-D  Fourier  transform 
of  the  modified  measured  data  to  the  3-D  Fourier  transform  of  a(r ) 
evaluated  over  a  semi-ellipsoid  oriented  at  angle  0.  As  the  measure¬ 
ment  angle  6  is  varied  from  0  to  27r,  two  overlapping  coverages  of  the 
Fourier  space  are  generated,  with  each  coverage  producing  one  of  the 
relationships  described  by  (20)  or  (21). 

Equations  (20)  and  (21)  yield  two  identical  values  of  Fkii'a^  ^'z) 
when  the  measured  data  are  consistent.  However,  when  the  measured 
data  contain  noise,  (20)  and  (2 1 )  will  generally  produce  different  values 
of  These  two  values  can  be  combined  to  obtain  a  final 

estimate  of  Pfc(i^Q,  J^z)  that  has  a  reduced  noise  level  as 

+  (1  —  a)*(l/m,  *'2)]  M*:(  — (22) 

where  ujkii'm,  I'z)  is  a  generally  complex -valued  combination  coef¬ 
ficient.  The  superscript  “a)’*  indicates  that  I'i)  is  obtained 

by  use  of  a  combination  coefficient  ujkii'm,  I'z)^  If  Mjt(i^m.  I'z)  and 
jy.)  are  interpreted  as  a  random  variables,  then  for  a  given 
i^kii'm,  i^zjf  (22)  can  be  interpreted  as  an  estimation  method  for  ob¬ 
taining  the  ideal  sinogram.  Because  I'z)  may  be  any  com¬ 

plex-valued  function  of  Um,  I'z  and  fc,  (22),  in  effect,  represents  an  in¬ 
finite  class  of  estimation  methods.  An  estimate  of  p(C,  0o)  can  be 
obtained  by  taking  the  3-D  inverse  Fourier  transform  of  P^"^^  {Pa .  Pz). 
For  a  fixed  value  of  the  filtered  backprojection  (FBP)  algorithm  of 
X-ray  CT  can  be  employed  to  reconstruct  the  corresponding  transverse 
slice  of  «(r,  9,  z)  from  p(^,  r,  (f>o).  We  refer  to  the  combination  of 
(22)  to  estimate  pI"'^  (i^m,  Pz)  coupled  with  the  2-D  FBP  algorithm  to 
reconstruct  transverse  slices  of  a(r,  9,  z)  as  a  spherical-wave  E-C  re¬ 
construction  algorithm  for  3-D  DT. 

IV.  Numerical  Results 


The  four  roots  of  (16)  are  given  by 


(17) 

(18) 


Leti/„.  =  =  1,2,3, and4.Fort/^+i^l^  <  i/;;(l  +  l/x'^),the 

two  roots  Prn3  and  p'm4  are  complex -valued  and,  therefore,  not  physi¬ 
cally  meaningful.  A  plot  of  (17)  is  shown  in  Fig.  2.  In  the  plane-wave 
case  (x  =  1),  there  are  only  two  roots  that  are  given  by  [4] 


where  p'i  p^  <  2po. 

Therefore,  for  a  given  pair  of  Pa  and  pi  satisfying  0  <  pI  pi  < 
Po{l-\-l/x^),Pk{Pa,  i^z)  can  be  obtained  from  M A:  (i^m,  i^z)ati^mi  = 
Pm  and  Pz  as 

Pt(r„,  u',)=ll{uLu  «':)l*M^(./,nI,  I',) 

=h(t/:„,r;)]*^Mi.(«/„„.^.-)  (20) 


We  performed  numerical  simulations  to  demonstrate  the  spher¬ 
ical-wave  E-C  reconstruction  algorithm.  We  considered  a  mathematical 
phantom  comprised  of  two  different  (uniform)  spheres  whose  3-D 
Fourier  transforms  were  approximately  bandlimited  to  a  sphere  of 
radius  y/2Po .  Our  intention  was  not  to  test  the  validity  of  the  weak 
scattering  (Bom)  condition  (see  the  Discussion  Section),  but  rather 
to  demonstrate  that  the  spherical-wave  E-C  reconstruction  algorithm 
can  accurately  reconstruct  the  scattering  object  function  from  weakly 
scattered  data.  We,  therefore,  employed  (3)  and  (6),  along  with  the 
analytic  expression  for  the  3-D  Fourier  transform  of  the  spheres,  to 
calculate  noiseless  samples  of  over  a  128  x  128  detector 

array  at  128  view  angles  that  were  evenly  spaced  over  360°.  In 
generating  the  simulation  data,  we  assumed  a  scanning  geometry 
with  5  =  100  (arbitrary  units)  and  D  =  104.0816  that,  according 
to  (4),  yields  x  —  0.7.  In  order  to  simulate  the  stochastic  nature 
of  noisy  scattered  data,  we  created  a  second  data  set,  Usi^^z^d), 
where  the  measured  scattered  data  were  treated  as  samples  of  an 
uncorrelated  bivariate  Gaussian  stochastic  process.  When  generating 
Us{^,z,<i>),  the  mean  and  variance  parameters  describing  the  real 
(imaginary)  component  of  the  stochastic  process  were  set  equal  to  the 
real  (imaginary)  values  of  the  noiseless  data  Us{^yZ,4>),  Therefore, 
at  a  given  position  (^,  z,0)  in  the  data  space,  the  magnitude  of  noise 
contained  in  the  real  and  imaginary  components  of  Ui,(i,z,<l>)  was 
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(c) 


Fig.  3.  (a)  True  slices  through  the  numerical  phantom.  Transverse  slices 
reconstructed  from  (b)  noiseless  and  (c)  noisy  data  using  a  spherical-wave  E-C 
reconstruction  algorithm. 

proportional  to  the  magnitude  of  the  real  and  imaginary  components 
of  2,  p),  respectively. 

We  reconstructed  three  transverse  slices  of  the  scattering  object  by 
use  of  the  spherical-wave  E-C  algorithm  specified  by  Ukii'm.  = 
1  /2  in  (22).  It  can  readily  be  shown  that  when  the  data  noise  is  uncor¬ 
related,  ic?  =  1/2  is  a  statistically  optimal  choice  in  the  sense  that  it 
minimizes  the  variance  of  the  estimate  ) ,  which  in  turn,  results 

in  a  minimization  of  the  global  image  variance^  [4].  For  more  general 
noise  models,  it  is,  in  principle,  possible  to  derive  other  optimal  forms 
for  i/r).  The  true  images  of  the  chosen  transverse  slices  are 

shown  in  Fig.  3(a).  The  images  reconstructed  from  the  noiseless  data, 
shown  in  Fig.  3(b),  do  not  contain  any  artifacts  and  accurately  repre¬ 
sent  the  corresponding  true  slices  shown  in  Fig.  3(a).  The  same  images 
reconstructed  from  the  noisy  data  set  are  shown  in  Fig.  3(c).  Out  of  cu¬ 
riosity,  we  also  reconstructed  the  same  transverse  slices  by  use  of  our 
previously  developed  plane-wave  E-C  reconstruction  algorithm  (that 
assumes  the  measurement  geometry  corresponds  to  x  =  1)  and  the 
noiseless  data  set.  The  reconstructed  images,  shown  in  Fig.  4,  clearly 
contain  artifacts  and  distortions.  This  numerically  demonstrates  the  im¬ 
portance  of  properly  accounting  for  the  wavefield  curvature  in  3-D  DT. 

V.  DISCUSSION 

Previously,  we  developed  a  novel  class  of  reconstruction  algorithms 
for  3-D  DT  using  plane-wave  sources.  These  algorithms,  referred  to  as 
plane-wave  E-C  reconstruction  algorithms,  had  a  significant  computa¬ 
tional  advantage  over  the  conventional  3-D  FBPP  algorithm,  and  unlike 
the  3-D  DF  method,  did  not  require  an  explicit  3-D  interpolation  in  the 
Fourier  space  of  the  scattering  object. 

^The  global  image  variance  is  equal  to  the  local  variance  of  an  image  at  a 
point  summed  over  all  points  in  image  space. 


Fig,  4.  Transverse  slices  reconstructed  from  the  spherical-wave  data  fiinction 
(x  =  0.7)  using  a  plane-wave  DT  reconstruction  algorithm.  As  expected,  the 
images  contain  distortions. 


The  use  of  a  plane-wave  source  may  not  be  feasible  in  many  ex¬ 
perimental  situations,  and  it  may  be  more  convenient  to  interrogate  the 
scattering  object  using  a  diverging  spherical  wave  that  is  produced  by  a 
point  source.  In  this  paper,  we  have  developed  a  class  of  spherical-wave 
E-C  reconstruction  algorithms  for  DT  using  spherical -wave  sources 
and  measurement  geometries  that  satisly  the  paraxial  approximation 
[7].  The  spherical-wave  E-C  reconstruction  algorithms  can  be  viewed 
as  generalizations  of  the  plane-wave  E-C  reconstruction  algorithms, 
and  reduce  to  the  plane- wave  E-C  algorithms  in  the  special  case  x  =  1  • 
The  spherical-wave  E-C  reconstruction  algorithms  possess  the  same 
advantages  as  their  plane-wave  counterparts.  For  example,  to  recon¬ 
struct  an  image  volume  the  spherical-wave  E-C  algorithm  requires 
log  N  numerical  operations  while  the  spherical-wave  FBPP  al¬ 
gorithm  described  by  (7)  would  require  log  A'  numerical  opera¬ 
tions.  Unlike  DF  methods,  the  spherical-wave  E-C  algorithms  do  not 
require  an  explicit  3-D  interpolation  in  the  Fourier  space  of  the  scat¬ 
tering  object. 

The  spherical-wave  E-C  reconstruction  algorithms  have  been  de¬ 
veloped  using  the  first-order  Bom  (or  Rytov)  weak  scattering  approx¬ 
imation.  In  certain  applications,  the  weak  scattering  approximation 
may  not  be  valid  and  the  reconstructed  image  may  contain  artifacts. 
However,  the  spherical-wave  E-C  algorithms  provide  a  natural  frame¬ 
work  for  the  incorporation  of  higher-order  scattering  perturbation 
approximations  [11]-[13]  into  the  algorithms. 
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A  Comparison  of  Algorithms  for  Detection  of  Spikes  in 
the  Electroencephalogram 
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Markad  V.  Kamath* 


Identification  of  the  short  transient  waveform,  called  a  spike, 
in  the  cortical  electroencephalogram  (EEC)  plays  an  important  role  during 
diagnosis  of  neurological  disorders  such  as  epilepsy.  It  has  been  suggested 
that  artificial  neural  networks  (ANN)  can  be  employed  for  spike  detec¬ 
tion  in  the  EEC,  if  suitable  features  are  provided  as  input  to  an  ANN.  In 
this  paper,  we  explore  the  performance  of  neural  network-based  classifiers 
using  features  selected  by  algorithms  suggested  by  four  previous  investiga¬ 
tors.  Of  these,  three  algorithms  model  the  spike  by  mathematical  parame¬ 
ters  and  use  them  as  features  for  classification  while  the  fourth  algorithm 
uses  raw  EEC  to  train  the  classifier.  The  objective  of  this  paper  is  to  ex¬ 
amine  if  there  is  any  inherent  advantage  to  any  particular  set  of  features, 
subject  to  the  condition  that  the  same  data  are  used  for  all  feature  selection 
algorithms.  Our  results  suggest  that  artificial  neural  networks  trained  with 
features  selected  using  any  one  of  the  above  three  algorithms  as  well  as  raw 
EEC  directly  fed  to  the  ANN  will  yield  similar  results. 

Index  Terms— Classification,  EEC,  neural  networks,  spike  detection. 


methods  that  mimic  the  mental  process  of  the  neurologist  in  identifying 
spikes  [5],  [6],  [17].  Research  involving  a  variety  of  signal  processing 
algorithms  has  led  to  automated  systems  for  analyzing  and  locating 
spikes  in  EEG  recordings  lasting  several  hours  [4],  [14].  While  it  is  de¬ 
sirable  to  have  detection  of  spikes  and  transient  waveforms  in  contin¬ 
uous  on-line  EEG,  efficient  algorithms  for  accurate  off-line  detection 
of  spikes  and  transient  waveforms  have  been  studied  for  some  time  [5], 

[11] ,  [14],  [19],  [21]. 

Previous  implementations  have  approached  the  problem  using 
back-propagation  [3],  [1 1],  [21]  and  Kohonen’s  self-organizing  maps 

[12] .  We  examined  four  methods  that  use  features  derived  from 
wavelet  transforms  (WTs)  [9],  autoregressive  (AR)  modeling  [18], 
context  parameters  [20]  and  an  ANN  trained  using  raw  EEG  [15]. 
These  algorithms  were  chosen  specifically  for  this  study  due  to  their 
high  spike  detection  rates  and  their  ability  to  reject  false  alarms. 

II.  Methods 

A.  Definition  of  Spikes 

Epileptic  EEG  contains  transient  waveforms  called  spikes  (or  spike 
discharge)  and  includes  the  following  types  of  waveforms  (Fig.  1). 

—  Spike:  lasts  20-70  ms. 

_ Sharp  wave:  lasts  70-200  ms  although  not  as  shaiply  contoured 

as  a  spike. 

—  Spike  and  wave  complex  :  A  spike  is  followed  by  a  slow  wave. 
If  they  occur  at  rates  below  3  Hz,  they  are  called  spike-and-slow 
wave  complexes. 

_ Polyspikes:  Multiple  spike  complexes;  several  spikes  occur  in 

sequence. 

_ Polyspike — and — slow  wave:  polyspikes  followed  by  a  slow 

wave. 

Fig.  1  shows  examples  of  a  spike,  spike  and  a  slow  wave,  polyspikes 
in  addition  to  normal  EEG,  and  EEG  contaminated  by  eye  blinks  and 
muscle  artifacts. 


I.  Introduction 

Epilepsy  is  characterized  by  sudden  recurrent  and  transient  distur¬ 
bances  of  mental  function  and/or  movements  of  the  body  that  result 
from  excessive  discharging  of  groups  of  brain  cells  [10],  [13].  Patients 
who  are  suspected  of  having  epileptogenic  foci  in  their  brain  are 
subjected  to  an  electroencephalography  (EEG)  recording  in  the  neu¬ 
rophysiology  laboratory.  In  clinical  neurological  practice,  detection 
of  abnormal  EEG  activity  plays  an  important  role  in  the  diagnosis  of 
epilepsy  [10].  It  is  generally  accepted  that  spikes  (often  called  “spike 
discharges”),  a  kind  of  transient  waveform(s)  present  in  human  EEG, 
have  a  high  correlation  with  seizure  occurrence.  Therefore,  detection 
of  spikes  in  the  EEG  plays  a  key  role  in  the  diagnosis  of  the  disease. 

Over  the  past  several  decades,  physicians  have  developed  empir¬ 
ical  techniques  which  help  them  identify  episodes  of  abnormal  signal 
components,  including  the  spike  discharges.  Such  expertise  has  led  to 


Manuscript  received  July  23,  200 1 ;  revised  November  22,  2002.  This  work 
was  supported  in  part  by  the  deGroote  Foundation  and  in  part  by  Natural  Sci¬ 
ences  and  Engineering  Research  Council  of  Canada  (NSERC).  Asterisk  indi¬ 
cates  corresponding  author. 

C.  C.  C.  Pang  is  with  the  Department  of  Electrical  and  Computer  Engineering, 
McMaster  University,  Hamilton,  ON  L8N  3Z5,  Canada. 

A.  R.  M.  Upton  and  G.  Shine  are  with  the  Department  of  Medicine.  McMaster 
University,  Hamilton,  ON  L8N  3Z5,  Canada. 

♦M.  V.  Kamath  is  with  the  Department  of  Medicine  and  the  Department 
of  Electrical  and  Computer  Engineering  and  Medicine,  McMaster  University, 
Room  3E25,  Health  Sciences  Bldg.,  Hamilton,  ON  L8N  3Z5,  Canada  (e-mail: 
kamathm@mcmail.cis.mcmaster.ca). 

Digital  Object  Identifier  10.1 1097rBME.2003. 809479 


B.  Feature  Extraction  Algorithms 

MathematicaT  modeling  of  spikes  has  enabled  researchers  a  wide 
range  of  methods  to  characterize  a  spike.  We  present  four  principal 
procedures  developed  by  previous  researchers  to  achieve  this  purpose. 

I)  Tarassenko's  Algorithm:  Tarassenko  et  al  [18]  considered 
both  time-domain  parameters  and  frequency-domain  parameters  for 
the  characterization  of  the  EEG  signal.  The  proposed  time-domain 
parameters  are  as  follows  (Fig.  2): 
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where,  Xav  is  the  average  amplitude  of  the  signal  within  an  epoch. 

The  latter  two  parameters,  which  were  proposed  by  Hjorth  [8]  and 
Walmsley  [19],  can  be  used  to  characterize  the  slope  and  slope  spread 
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ABSTRACT:  It  is  widely  believed  that  measurements  from  a  full 
angular  range  of  27r  are  generally  required  to  exactly  reconstruct  a 
complex-valued  refractive  index  distribution  in  diffraction  tomography 
(DTI.  in  this  work,  we  developed  a  new  class  of  minimal-scan  recon¬ 
struction  algorithms  for  DT  that  utilizes  measurements  only  over  the 
angular  range  0  <  ^  37r/2  to  perform  an  exact  reconstruction. 

These  algorithms,  referred  to  as  minimal-scan  estimate-combination 
(MS-E-C)  reconstruction  algorithms,  effectively  operate  by  transform¬ 
ing  the  DT  reconstruction  problem  into  a  conventional  x-ray  CT  re¬ 
construction  problem  that  requires  inversion  of  the  Radon  transform. 
We  performed  computer  simulations  to  compare  the  noise  and  nu¬ 
merical  properties  of  the  MS-E-C  algorithms  against  existing  filtered 
backpropagation-based  algorithms.  ©  2002  Wiley  Periodicals,  Inc.  Int  J 
Imaging  SystTechnol.  12,  84-91,  2002;  Published  online  In  Wiley  InterScience 
(www.Interscience. wiley.com).  DOl  10.1002/ima.10014 

Key  words:  topographic  reconstruction;  diffraction  tomography; 
wavefield  inversion 


I.  INTRODUCTION 

In  diffraction  tomography  (DT),  a  scattering  object  is  interrogated 
using  a  diffracting  acoustical  or  electromagnetic  wavefield,  and  the 
scattered  wavefield  around  the  object  is  measured  and  used  to 
reconstruct  the  refractive  index  distribution  of  the  scattering  object. 
There  are  numerous  potential  applications  of  DT  that  can  be  found 
in  various  scientific  fields  (Andre  et  al.,  1995;  Tabbara  et  al.,  1988; 
Mueller  et  al.,  1979;  Kino,  1979;  Devaney,  1984;  Robinson,  1984). 
Recently,  there  has  also  been  considerable  interest  in  using  DT  to 
perform  coherent  x-ray  imaging  using  third-generation  synchrotron 
sources  (Cheng  and  Han,  2001).  Unlike  the  x-rays  used  in  computed 
tomography  (CT)  that  travel  along  straight  lines,  the  radiation  em¬ 
ployed  in  DT  has  to  be  treated  in  terms  of  wavefronts  and  fields 
scattered  by  inhomogeneities  in  the  object.  In  DT,  the  interaction 
between  the  incident  wavefield  and  the  object  medium  is  governed 
by  the  inhomogeneous  Helmholtz  equation.  Using  a  weak-scattering 
approximation,  the  inhomogeneous  equation  can  be  analytically 
solved  (Wolf,  1969;  Mueller  et  al.,  1979)  to  obtain  a  linear  relation¬ 
ship  between  the  scattered  field  and  the  refractive  index  distribution. 


Correspondence  lo:  Mark  A.  Anastasio,  Ph.D..  Pritzker  Institute  of  Medical  Engi¬ 
neering,  Illinois  Institute  of  Technology,  10  We.st  32nd  Street.  Chicago.  IL  60616-3793. 
(V)  312-567-.3926,  (F)  312-.S67-5707,  Email:  anastasio® iit.cdu. 


This  relationship  has  been  used  to  develop  DT  reconstruction  algo¬ 
rithms  such  as  the  well-known  filtered  backpropagation  (FBPP) 
algorithm  (Devaney,  1982),  which  is  a  generalization  of  the  filtered 
backprojection  (FBPJ)  algorithm  of  x-ray  CT. 

It  is  widely  believed  that  measurements  from  a  full -angular 
range  of  Itt  around  the  scattering  object  are  generally  required  to 
exactly  reconstruct  a  complex-valued  refractive  index  distribu¬ 
tion  (Devaney,  1982).  However,  we  have  recently  revealed  that 
one  needs  measurements  only  over  the  angular  range  0  ^  ^ 

37r/2  to  perform  an  exact  reconstruction,  and  we  developed 
minimal -scan  filtered  backpropagation  (MS-FBPP)  algorithms  to 
achieve  this  (Pan  and  Anastasio,  1999).  A  useful  characteristic  of 
the  MS-FBPP  algorithms  is  their  ability  to  decrease  the  data 
acquisition  time  by  at  least  25%  over  conventional  (full-scan) 
algorithms.  They  can  also  reduce  artifacts  due  to  movement  in  or 
temporal  fluctuations  of  the  scattering  object.  Furthermore,  in 
certain  practical  situations,  it  may  be  impossible  to  acquire 
measurements  over  a  full  2'it  angular  range. 

A  new  class  of  reconstruction  algorithms  has  recently  been 
developed  for  full-scan  DT  (in  other  words,  DT  employing  mea¬ 
surements  over  a  full  27r  angular  range).  These  algorithms,  referred 
to  as  estimate-combination  (E-C)  reconstruction  algorithms  (Pan, 
1998;  Anastasio  and  Pan,  2000b;  Anastasio  and  Pan,  2000a),  effec¬ 
tively  operate  by  transforming  the  DT  reconstruction  problem  into  a 
conventional  x-ray  CT  reconstruction  problem  that  can  be  efficiently 
solved  using  the  filtered  backprojection  (FBPJ)  algorithm.  The  E-C 
reconstruction  algorithms  are  more  computationally  efficient  than 
the  FBPP  algorithm,  and  also  provide  a  flexible  framework  for 
imposing  unbiased  regularization. 

Because  the  E-C  reconstruction  algorithms  involve  a  Fourier 
series  expansion  of  the  data  function  that  is  acquired  over  the 
angular  range  0  <  <j>  ^  27r,  they  cannot  be  applied  directly  to  the 
minimal-scan  problem  where  measurements  are  only  acquired  over 
the  angular  range  0  <  (f>  ^  37r/2.  Because  of  the  potential  advan¬ 
tages  of  the  E-C  reconstruction  algorithms,  it  is  important  to  gen¬ 
eralize  them  to  the  minimal-scan  situation.  In  this  work,  we  devel¬ 
oped  minimal-scan  E-C  (MS-E-C)  reconstruction  algorithms  for  DT. 
We  performed  computer  simulations  to  compare  the  noise  and 
numerical  properties  of  the  MS-E-C  and  MS-FBPP  algorithms.  Our 
results  quantitatively  demonstrate  that  the  MS-E-C  algorithms  pos- 
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Scattered  data 


Figure  1 .  The  classical  scanning  geometry  of  2D  DT.  The  insonifylng 
plane  wave  propagates  along  the  iq  axis,  and  the  scattered  wave  field 
is  measured  along  the  line  tj  I.  Full-scan  and  minimal-scan  data  sets 
are  obtained  by  varying  the  measurement  angle  <f)  between  0  and  27r 
or  between  0  and  877/2,  respectively. 

sess  statistical  and  numerical  properties  superior  to  those  of  the 
MS-FBPP  algorithms. 

II.  BACKGROUND 

A.  The  Fourier  Diffraction  Projection  Theorem.  In  two-di¬ 
mensional  (2D)  DT  employing  the  classical  scanning  configura¬ 
tion,  as  shown  in  Figure  1,  the  scattering  object  is  illuminated  by 
monochromatic  plane-wave  radiation  of  frequency  and  the 
transmitted  wavefield  is  measured  along  the  f  axis  oriented  at  a 
mea.surement  angle  (J>,  at  a  distance  rj  =  I  from  the  origin.  From 
measurements  of  the  scattered  wavefield  obtained  at  various 
angles  </>,  one  seeks  to  reconstruct  the  scattering  object  function 
a(r),  which  is  related  to  the  refractive  index  distribution  n(r,  6)  by 
a(r.  0)  =  /r(r,  0)  -  1. 

At  a  measurement  angle  cl>.  the  scattered  data  are  measured  along 
the  line  tj  =  /,  as  shown  in  figure  1.  Let  </>)  to  denote  the 

ID  Fourier  transform  of  the  measured  scattered  data  with  respect  to 
For  convenience,  we  define  a  modified  ID  Fourier  transform  of 
the  scattered  data  as 

=  (1) 

where  v'  =  Vvf,  -  vi  and  |v„J  s  Vo.  The  quantities  on  the 
right-hand  side  of  equation  1  are  known  or  can  be  measured. 
Therefore,  we  will  treat  M{v„,  </>)  as  a  measurable  data  function. 
Under  the  Bom  approximation  (Mueller  et  al.,  1979),  the  Fourier 
diffraction  projection  (FDP)  theorem  (Mueller  et  al.,  1979)  can  be 
derived,  which  is  mathematically  stated  as 

4.)  =  J  I 

=  0  if|vj>^), 

(2) 


where  the  polar  coordinates  (r,  0)  and  the  rotated  coordinates  (^,  tj) 
are  related  through  ^  =  r  cos{4>  —  6)  and  tj  =  —  r  sin(^  ~  0).  The 
FDP  theorem  indicates  that  (f))  provides  the  values  of  the  2D 

Fourier  transform  of  air)  along  the  semi-circular  arc  AOB  of  radius 
i/q,  as  shown  in  figure  2. 

B.  Minimal-Scan  Filtered  Backpropagation  Algorithms. 

The  widely  used  filtered  backpropagation  (FBPP)  algorithm  (Dev- 
aney,  1982)  is  mathematically  expressed  as 

aCr)  =  ^  7  K\M(v„,  dv^  d4>. 

<lf=0  ^  Vm=-tin 

(3) 

where  =  jiV^Vo  “  ^ben  Vq  — >  the  FBPP 

algorithm  reduces  to  the  filtered  backprojection  (FBPJ)  algorithm  of 
x-ray  CT.  The  FBPP  algorithm  generally  requires  full  knowledge  of 
(t>)  in  the  data  space  ‘W*  =  <  Vq,  0  ^  d)  ^  2tt\,  for 

exact  reconstruction  of  the  generally  complex-valued  object  func¬ 
tion.  We  will  refer  to  such  full  knowledge  of  A/(  <!>)  as  a  full-scan 

data  set.  The  fiill-scan  data  space  can  be  decomposed  into  four 
subspaces,  si,  95,  and  Q),  where  si  =  [|v„|  <  Vq,  ,  0  <  <l)  ^ 

2a  +  7r/2],  95  =  <  i/q,  7r/2  +  2a  <  d)  ^  tt  +  2aj,  ^ 

=  :<  i/q,  TT  +  2a  <  37r/2]  and  Vo> 

37r/2  <  d>  2.7r],  where  a  =  sgn(y„,)arcsinVr^^^^^/(2rJ.  A 
schematic  of  this  partitioning  of  the  data  space  is  given  in  figure  3. 
Using  the  FDP  theorem,  it  can  be  shown  (Pan  and  Kak,  1983)  that 
Miv„,,  d>)  -  d>  +  "  2a).  Therefore,  the  information 

contained  in  subspace  si  is  redundant  to  that  contained  in  subspace 
and  the  information  contained  in  subspace  is  redundant  to  that 
contained  in  subspace  Q).  We  have  demonstrated  (Pan  and  Anasta- 
sio,  1999)  that  it  is  possible  to  exactly  reconstruct  the  object  function 
using  only  knowledge  of  M{  v„,,  d>)  in  the  subspace  jtf  =  U 


Figure  2.  The  FDP  theorem  states  that  d>)  is  equal  to  the  2D 
Fourier  transform  of  a(f)  along  a  semi-circle  AOB  that  has  a  radius  of 
1^0  and  is  centered  at  =  vq. 
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Figure  3.  The  subspaces  ^  and  2)  in  the  complete  DT  data 
space. 


y  eg  =  [|^^|  <  i/y,  0  :<  r<  37r/2],  which  we  refer  to  as  a  a 
minimal-scan  data  set. 

Although  the  minimal-scan  data  set  il  contains  all  of  the  infor¬ 
mation  necessary  for  exact  reconstruction  of  the  scattering  object 
function,  the  redundant  information  contained  in  the  subspaces  si 
and  ^  needs  to  be  properly  normalized  in  the  reconstruction  process 
(Pan  and  Anastasio.  1999).  The  MS-FBPP  algorithms  operate  by 
first  normalizing  such  partially  redundant  information  by  generating 
an  appropriately  weighted  minimal-scan  data  set  Af <!))  and 
subsequently  using  the  FBPP  algorithm  described  by  equation  3 
(scaled  by  a  factor  of  2),  to  exactly  reconstruct  the  image.  The 
weighted  minimal-scan  data  set  is  given  by  (Pan  and  Anastasio. 
1999) 


<#>)  =  (4) 

where  4>)  is  a  function  of  and  <f>,  which  satisfies 

wi  v,„,  (f))  +  w{  —  v,„,  </)  4-  TT  —  2a)  =  1  (5a) 

in  complete  data  space  TV", 


<t))  =  1  (5b) 

in  subspace  2ft.  and 

(/>)  =  0  (5c) 


in  subspace  05.  Although  the  forms  of  <j>)  in  subspaces  2ft  and 
9)  are  completely  specified  by  equations  5b  and  5c,  respectively,  the 
explicit  forms  of  w(  (t>)  in  subspaces  .<4  and  ^  are  unspecified  for 
the  moment.  In  principle,  one  can  choose  different  >v(v„,  <^))  in 
subspaces  ^4  and  ^  as  long  as  these  w(r,„,  <^>)  satisfy  equation  5a. 

III.  MINIMAL-SCAN  ESTIMATE-COMBINATION 
ALGORITHMS 

The  previously  derived  (full-scan)  E-C  reconstruction  algorithms  are 
more  computationally  efficient  than  the  (full-scan)  FBPP  reconstruc¬ 
tion  algorithms,  which  involve  a  depth-dependent  filtering  operation 
(backpropagation).  Accordingly,  we  expect  that  the  MS-FBPP  algo¬ 
rithms,  which  use  the  FBPP  algorithm  to  reconstruct  the  final  image 
from  the  weighted  data  function  <t>),  will  also  be  less 

computationally  efficient  than  the  E-C  reconstruction  algorithms. 
Because  they  will  involve  fewer  and  less  complicated  numerical 
operations,  we  also  expect  that  the  MS-E-C  algorithms  will  be 
minimize  the  propagation  of  data  noise  and  errors  as  compared  to 
the  MS-FBPP  algorithms.  The  full-scan  E-C  reconstruction  algo¬ 
rithms  (Pan,  1998;  Anastasio  and  Pan,  2000b)  involve  a  Fourier 
series  expansion  of  the  data  function  <t>).  which  requires 

knowledge  of  M(v„,  <j>)  over  an  angular  range  of  2^-,  and  therefore 
can  not  be  directly  applied  to  the  minimal-scan  data  set  containing 
only  measurements  in  the  range  0  ^  <j>^  317/3.  Below,  we  develop 
minimal-scan  E-C  (MS-E-C)  reconstruction  algorithms  that  can  be 
directly  applied  to  the  minimal-scan  data  set. 

A.  The  Radon  Transform.  Let  pit  </>)  and  P^iK)  denote  the 
Radon  transform  of  «(r,  6)  and  its  2D  Fourier  transform,  respec¬ 
tively.  (Here,  the  2D  Fourier  transform  is  actually  a  ID  Fourier 
transform  with  respect  to  ^  and  a  ID  Fourier  series  with  respect  to 
(j).)  From  knowledge  of  pit  <^>)^  or  equivalently,  P*(v„),  one  can 
reconstruct  air,  0)  by  use  of  the  computationally  efficient  and 
numerically  stable  FBPl  algorithm,  which  is  given  by 

a{r,  0)  =  i  r  I"  i 

(6a) 

For  theoretical  convenience,  the  FBPJ  algorithm  can  also  be  ex- 
pressed  as 


a(r,  e)  =  Ztt  f  \  Fi(i'„)^‘''A(27rv„r)v„  (6b) 

where  Ji,i  •  )  is  a  Bessel  function  of  the  first  kind.  The  MS-E-C 
algorithms  will  operate  by  estimating  or  equivalently  pit 

if)),  from  the  minimal-scan  data  set,  and  using  the  FBPJ  algorithm  to 
reconstruct  the  final  image  air,  B). 

B.  Derivation  of  the  MS-E-C  Algorithms.  Consider  a  given 
weighting  function  wiv^,  <#>)  in  equation  4.  The  corresponding 
MS-FBPP  algorithm  can  be  expressed  as  (Pan  and  Anastasio,  1999) 


a{r,  0)  =  2 


^  I  dv„  d<l>. 


(7) 
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Let  M\(v„)  denote  the  Fourier  series  expansion  of  4>).  One 

can  re-express  equation  7  as 

air.  6) 


=  -  I"  I  E  Mlivje^-*  di’„d4>. 

(8) 

Using  the  definition  y{v„)  =  and  separating  the  contribution  to 
the  integral  from  positive  and  negative  equation  8  can  be 
re-written  as 


a{r,  e)  = 


i 


I/' 


« 

E 


Ml(v„)  dv,„  d4> 


X  E  M[(vjdv„dd>.  (9) 


Changing  v„  to  -  v,„  in  the  second  term  in  equation  9  and  grouping 
d)-dependent  terms  yields 


fi(r,  0)  =  2 


d^ 


X  2 


2  M[(-v„)dv^, 


(10) 


and  therefore  the  Radon  transform  of  the  scattering  object  function 
a(r,  0)  can  be  estimated  from  the  appropriately  weighted  minimal- 
scan  data  set.  The  use  of  equation  12  to  estimate  P  ifvf)  coupled  with 
the  2D  FBPJ  algorithm  to  reconstruct  a{r,  6)  is  referred  to  as  a 
MS-E-C  reconstruction  algorithm.  In  practice,  the  FBPJ  algorithm 
described  by  equation  6a  requires  knowledge  of  P*(  vj  for  evenly 
spaced  values  of  spanning  the  range  -  V2  ^  V2  Vq 

in  order  to  be  efficiently  implemented  using  the  fast  Fourier  trans¬ 
form  (FFT).  In  this  case,  the  consistency  condition  (Deans,  1983) 
p^{vj  =  (-‘l)*P*(“J^ci)  can  be  employed  to  obtain  values  of 
Pki^a)  for  negative  v^. 

IV.  NUMERICAL  SIMULATIONS 

We  performed  simulation  studies  to  evaluate  and  compare  the  nu¬ 
merical  and  statistical  properties  of  images  obtained  by  use  of  the 
MS-E-C  and  MS-FBPP  reconstruction  algorithms. 

A.  Measurement  Data.  We  investigated  the  statistical  proper¬ 
ties  of  the  reconstruction  algorithms  under  near-ideal  conditions  by 
employing  a  single  component  scattering  object  that  exactly  satisfied 
the  (first-order)  Bom  approximation.  The  propagation  of  determin¬ 
istic  artifacts  by  the  reconstruction  algorithms  under  less-than-ideal 
conditions  was  investigated  by  employing  a  two  component  .scatter¬ 
ing  object  that  introduced  strong  and  multiple-scattering  effects  into 
the  measurement  data. 

A.  7  Single-Scattered  (Bom)  Data  and  Noise  Model.  The  scat¬ 
tering  object  function,  shown  in  figure  4,  was  taken  to  be  a  lossless, 
uniform  cylindrical  disk  with  a  diameter  of  30  pixels  that  was 
convolved  with  a  symmetric  Gaussian  function  with  a  standard 
deviation  of  0.2  pixel.  The  Fourier  transform  of  the  object  function 
was  therefore  approximately  bandlimited  to  a  circular  disk  of  radius 
V2  1^0  in  its  2D  Fourier  space.  It  was  assumed  that  the  scatterer  was 
weakly  scattering,  so  that  the  Born  approximation  may  reasonably 
be  taken  to  hold.  Therefore,  these  numerical  simulations  were  de¬ 
signed  to  investigate  the  statistical  properties  of  the  reconstruction 
algorithms  rather  than  the  weak  scattering  model.  Minimal-scan  data 
sets  were  generated  by  using  the  FDP  theorem  to  calculate  simulated 
scattered  field  data  for  96  measurement  angles  </>  that  were  evenly 
spaced  between  0  and  37r/2.  At  each  measurement  angle.  128 
samples  were  calculated  with  a  sampling  increment  =  1/2 
where  Vq  is  the  frequency  of  the  incident  plane  wave. 


The  integrals  in  the  curly  braces  of  the  first  and  second  terms  of 
equation  10  can  be  evaluated  (Metz  and  Pan  1995)  yielding  the 
expressions  f,{2Trr'\/ vf„  ~  v^)  and 

27r(-  1  )y7(i^,„)"V^%(27rrVv^,  -  ^),  respectively.  Using 
this  result  and  the  change  of  variables  (which 

implies  v„dv^,  =  vjv  v„,dv^).  we  arrive  at 

"  r 

o(r,  6)  =  IT  E  /  WM'Av^)  +  (-l)‘y'X(-»'™)3 

X  (2Trv^r) v„  dv^,  (11) 

where  v,„  =  r^Vl  -  vIHvq.  Comparison  of  Eqns.  1 1  and  6b 
indicates  that,  for  0  ^  ^ 

P^M  =  J  [yA/KO  +  (-ify-'MK-vJl  (12) 


Figure  4.  The  original  scattering  object  function  was  formed  by 
convolving  a  uniform  circular  disk  with  a  diameter  of  30  pixels  with  a 
circularly  symmetric  Gaussian  function  with  a  standard  deviation  of 
0,2  pixels. 
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Figure  5.  Images  reconstructed  from  noiseless  data  using  the  (a) 
MS-E-C  and  (b)  MS-FBPP  reconstruction  algorithms. 


To  simulate  the  effects  of  data  noise,  we  treated  the  scattered  data 
<#>)  as  a  complex  stochastic  process  with  a  real  and  an  imag- 
inary  component  denoted  by,  <f>)  and  <#>),  respectively. 

Let  u*''-  =  and  where  and 

are  the  means  of  u‘/’  and  u''\  respectively.  The  statistics  of  the 
deviates  Au*''*  and  AuJ'*  are  described  by  the  circular  Gaussian 
model: 


p(Aii''\  Au'")  = 


1  r  1  /A«'^>'  A«"’'\' 


(13) 


where  a;  and  af  are  the  variances  of  (#>)  and  <<>), 

respectively. 

A.2  Multiple-Scattered  Data,  To  investigate  the  impact  of 
strong  and  multiple-scattering  effects  on  the  performance  of  the 
MS-E-C  and  MS-FBPP  algorithms,  which  implicitly  assume  weak 
scattering  conditions,  we  employed  the  two  component  scattering 
object  shown  in  figure  8.  This  scattering  object  was  compo.sed  of 
two  uniform  cylinders  with  radius  3A  whose  centers  were  separated 
by  7A.  (A  =  wavelength  of  incident  plane  wave.)  The  refractive 
index  values  of  the  cylinders  were  varied  as  described  below. 
Twersky’s  theory  of  multiple  scattering  (Ishimaru,  1978)  was  used 
to  calculate  measurement  data  that  contained  second-order  scatter¬ 
ing  effects.  The  first-order  contributions  to  the  measurement  data 
were  obtained  by  considering  the  interaction  of  the  incident  plane- 
wave  with  each  cylinder,  assuming  the  other  cylinder  to  be  absent. 
Note  that  this  was  an  exact  calculation  that  did  not  rely  on  the  Bom 
approximation.  The  second-order  contributions  to  the  measurement 
data  were  obtained  by  calculating  the  scattering  component  created 


(«)  W 

Figure  6.  Images  reconstructed  from  noisy  and  weakly  scattered 
data  using  the  (a)  MS-E-C  and  (b)  MS-FBPP  reconstruction  algo¬ 
rithms.  The  noisy  data  were  generated  using  =  u,  =  0.5  in  Eqn.  13. 


Figure  7.  A  plot  of  the  local  variance  obtained  from  the  MS-FBPP 
reconstructed  images  divided  by  the  local  variance  obtained  from  the 
MS-E-C  reconstructed  images. 


when  the  incident  wave  interacts  with  one  cylinder  and  subsequently 
scatters  off  of  the  second  cylinder  before  being  measured  (Azimi 
and  Kak,  1983).  Minimal-scan  data  sets  were  generated  containing 
96  measurement  angles  (j)  that  were  evenly  spaced  between  0  and 
3tt/2.  At  each  measurement  angle,  128  samples  were  calculated  with 
a  sampling  increment  —  1/2^^. 

B.  Implementation  Details.  Both  the  MS-E-C  and  MS-FBPP 
reconstruction  algorithms  require  the  computation  of  the  weighted 
minimal-scan  data  set  that  is  defined  in  equation  4.  This  was  ac¬ 
complished  by  using  the  weighting  function 


f 


w(v,„,  (#))  = 


1  TT 

:  ^  +  2a  ^  <f>  ^  TT  +  2a 

1  377 

-:7r  +  2a^<f)^'y 

377 

0  :  —  <  ^  ^  277, 


(14) 


Figure  8.  The  geometry  of  the  two  component  scattering  object.  We 
considered  the  cases  where  the  refractive  Index  of  the  uniform  cylin¬ 
ders  took  the  values  n(r)  =  1.01, 1.05,  and  1.08. 
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where  |i^„J  ^  Vq.  This  weighting  function  satisfies  the  requirements 
described  by  equation  5. 

MS-E-C  Algorithm.  From  the  uniformly  sampled  values  of  the 
scattered  field  hence  can  be  deter¬ 

mined  at  uniformly  spaced  values  of  v„f.  (This  calculation  can  be 
performed  using  the  FFT  algorithm.)  However,  the  uniformly 
spaced  values  of  at  which  is  known  do  not  generally 

correspond  to  the  uniformly  spaced  values  of  at 

which  we  are  required  to  evaluate  For  each  of  65  evenly 

spaced  values  of  spanning  the  range  0  <  ^  V2  Vq,  we  used 

linear  interpolation  to  determine  the  values  of  A/i  (v,„  = 
y/vl  +  vl)  and  A/i  ( =  -Vvl  -b  v^)  from  the  sampled  values 
of  and  which  were  subsequently  used  in  equa¬ 
tion  13  To  evaluate  Por  each  value  of  k.  zero-padding 

interpolation  was  employed  to  increase  the  sampling  density  along 
the  axis  of  by  a  factor  of  three  in  order  to  increase  the 

accuracy  of  the  interpolation  operation.  The  consistency  condition 
/5^(  pj  =  (-\  )*/>^(  -  vj  was  employed  to  obtain  the  64  samples  of 
Pf,{vJ  for  spanning  the  range  -V2  <  0.  In  our 

implementation  of  the  FBPJ  algorithm,  an  unapodized  ramp  filter 
was  used.  The  interpolation  necessary  to  align  the  backprojected 
data  onto  a  128X128  pixel  discrete  image  matrix  was  performed 
using  bilinear  interpolation. 

MS’FBPP  Algorithm.  In  the  MS-FBPP  algorithm,  at  each  mea¬ 
surement  angle  <f>.  M{  </>)  was  multiplied  by  the  depth-dependent 
filter  Vq/v'  |  for  each  of  128  discrete  values  of  tj.  For  each 

value  of  7),  the  filtered  data  was  zero-padded  to  a  length  of  182 
samples,  upon  which  the  inverse  FFT  was  employed  to  transform  the 
filtered  data  to  the  corresponding  depth  (value  of  rj)  in  image  space. 
This  ensured  that  the  pixel  size  of  the  images  reconstructed  using  the 
MS-E-C  and  MS-FBPP  algorithms  were  equivalent.  The  interpola¬ 
tion  necessary  to  align  the  backpropagated  data  onto  a  128x128 
pixel  discrete  image  matrix  was  performed  using  bilinear  interpola¬ 
tion. 

C.  Simulation  Studies 

C./  Single-Scattered  (Bom)  Data  Case.  To  study  the  noise 
properties  of  the  reconstructed  images  quantitatively,  we  utilized  the 
noiseless  data  set  described  in  Section  IV-A.l  along  with  the  noise 
model  in  equation  13  with  a,.  =  <t,  ==  0.5  to  generate  =  250 
noisy  minimal-scan  data  sets.  The  MS-E-C  and  MS-FBPP  algo¬ 
rithms  were  used  to  reconstruct  sets  of  250  noisy  images  from  these 
noisy  minimal-scan  data  sets.  The  local  image  variance  was  calcu¬ 
lated  empirically  from  the  N  sets  of  reconstructed  images  as 


(M 

Figure  10.  Images  of  the  two  component  scattering  object  corre¬ 
sponding  to  n(f)  =  1 .05,  reconstructed  by  use  of  the  (a)  MS-E-C  and 
(b)  MS-FBPP  reconstruction  algorithms. 


Var{a(r)}  =  a,(rf  -  ^  j  j’ 

where  a,(r)  is  the  ith  image  obtained  with  either  the  MS-E-C  or 
MS-FBPP  reconstruction  algorithm. 

C.2  Multiple-Scattering  Case.  To  assess  the  impact  of  multiple¬ 
scattering  effects  on  the  performance  of  the  MS-E-C  and  MS-FBPP 
algorithms,  we  employed  the  two  component  scattering  object  de¬ 
scribed  in  Section  IV-A.2.  We  considered  the  three  cases  where  the 
cylinders  had  refractive  index  values  of  n(r)  =  1.01,  1.05,  and 
1.08.  For  each  value  of  «(r),  we  used  the  MS-E-C  and  MS-FBPP 
algorithms  to  reconstruct  the  image  from  the  corresponding  data  set. 
The  percentage  of  cumulative  error  in  the  reconstructed  images  was 
quantified  using  the  metric 

|a(r)  -  dr 

Ero,  =  - X 

I  \ci,ru.(f)\^  dr 

ROl 

where  a,^^^(r)  is  the  true  scattering  object  function  and  the  subscript 
‘ROr  denotes  that  the  error  was  calculated  over  a  64X64  pixeP 
region  of  interest  containing  the  scattering  objects. 


Figure  9.  Images  of  the  two  component  scattering  object  corre¬ 
sponding  to  n(f)  =  1.01,  reconstructed  by  use  of  the  (a)  MS-E-C  and 
(b)  MS-FBPP  reconstruction  algorithms. 


(*)  (b) 

Figure  11.  Images  of  the  two  component  scattering  object  corre¬ 
sponding  to  /?(F)  =  1 .08,  reconstructed  by  use  of  the  (a)  MS-E-C  and 
(b)  MS-FBPP  reconstruction  algorithms. 
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Table  I.  Error  values  of  the  reconstructed  images  shown  in  Figs.  9-11. 


Contrast 

MS-E-C  Error 

MS-FBPP  Error 

1.01 

12.03 

13.68 

1.05 

51.26 

53.70 

1.08 

115.25 

119.13 

V.  RESULTS 

A.  Single-Scattered  (Born)  Data  Case.  We  first  used  the  MS- 
E-C  and  MS-FBPP  algorithms  to  reconstruct  the  scattering  object 
function  using  the  simulated  noiseless  minimal-scan  data  set.  The 
images  reconstructed  using  the  MS-E-C  and  MS-FBPP  algorithms 
are  displayed  in  figures  5a  and  5b,  respectively.  It  is  observed  that, 
in  the  absence  of  noise,  both  the  MS-E-C  and  MS-FBPP  algorithms 
can,  with  high  fidelity,  reconstruct  the  original  scattering  object 
function  from  the  minimal-scan  data  set. 

Using  one  of  the  noisy  minimal-scan  data  sets,  we  used  the 
MS-E-C  and  MS-FBPP  algorithms  to  reconstruct  the  scattering 
object  function.  The  noisy  images  reconstructed  using  the  MS-E-C 
and  MS-FBPP  algorithms  are  displayed  in  figures  6a  and  6b.  re¬ 
spectively.  The  image  reconstructed  using  the  MS-E-C  algorithm 
(Fig,  6a)  appears  less  affected  by  the  data  noise  and  more  closely 
resembles  the  original  object  than  does  the  image  reconstructed 
using  the  MS-FBPP  algorithm  (Fig.  6b).  The  local  image  variance, 
which  was  empirically  calculated  from  the  two  sets  of  250  noisy 
images  reconstructed  using  the  MS-E-C  and  MS-FBPP  algorithms, 
quantitatively  confirms  this  observation.  Figure  7  is  a  plot  of  the 
local  variance  obtained  from  the  MS-FBPP  reconstructed  images 
divided  by  the  local  variance  obtained  from  the  MS-E-C  recon¬ 
structed  images.  Clearly,  the  ratio  of  the  variances  is  everywhere 
greater  than  one,  and  near  the  corners  of  the  reconstructed  image  is 
as  great  as  ten.  This  quantitatively  demonstrates  that  the  MS-E-C 
reconstruction  algorithms  are  less  susceptible  to  the  effects  of  data 
noise  than  are  the  MS-FBPP  reconstruction  algorithms. 

B.  Multiple-Scattering  Case.  Using  the  MS-E-C  and  MS- 
FBPP  algorithms  we  reconstructed  the  two  component  scattering 
objects  shown  in  figures  9-1  i,  which  correspond  to  the  cases  where 
the  cylinders  had  refractive  index  values  of  n{r)  -  1.01,  1 .05,  and 
1.08,  respectively.  In  each  case,  the  image  reconstructed  by  use  of 
the  MS-E-C  algorithm  (Figs.  9a-l  la)  appears  to  contain  less  pro¬ 
nounced  artifacts  than  does  the  image  reconstructed  by  use  of  the 
MS-FBPP  algorithm  (Figs.  9b-llb).  This  observation  is  confirmed 
by  Table  I,  which  shows  that  for  each  value  of  ti{r)  the  MS-E-C 
algorithm  produced  images  that  have  lower  error  values  than  the 
corresponding  images  produced  by  the  MS-FBPP  algorithm.  This 
quantitatively  demonstrates  that  the  MS-E-C  algorithms  are  less 
susceptible  to  multiple-scattering  effects  and  other  deterministic 
inconsistencies  than  are  the  MS-FBPP  algorithms.  However,  as  one 
would  expect,  the  performance  of  both  algorithms  dramatically 
deteriorates  as  the  refractive  index  values  increases  and  the  Bom 
condition  (Chen  and  Stamnes,  1998)  is  severely  violated. 

VI.  DISCUSSION 

Previously  we  have  shown  (Pan  and  Anastasio,  1999)  that,  in  DT 
employing  the  2D  classical  scanning  geometry,  the  minimal-scan 
data  set  acquired  using  view  angles  only  in  [0,  37r/2J  contains  all  of 
the  information  necessary  to  exactly  reconstruct  the  scattering  object 
function.  We  subsequently  developed  a  class  of  MS-FBPP  algo¬ 


rithms  that  were  capable  of  exactly  reconstructing  the  scattering 
object  function  from  the  minimal-scan  data  set. 

In  this  work,  we  have  developed  a  novel  class  of  reconstruction 
algorithms  for  the  minimal-scan  DT  reconstruction  problem.  These 
algorithms,  referred  to  as  MS-E-C  reconstruction  algorithms,  have 
distinct  advantages  over  the  MS-FBPP  reconstruction  algorithms. 
Because  the  FBPJ  algorithm  used  by  the  MS-E-C  algorithms  does 
not  involve  a  depth-dependent  filtering,  the  MS-E-C  algorithms  are 
more  computationally  efficient  than  are  the  MS-FBPP  algorithms. 
More  importantly,  we  have  quantitatively  demonstrated  that  the 
MS-E-C  algorithms  are  less  susceptible  to  data  noise,  modeling 
errors  due  to  the  violation  of  weak  scattering  conditions,  and  other 
finite  sampling  effects  than  are  the  MS-FBPP  algorithms.  This  result 
is  consistent  with  the  observation  that  the  MS-FBPP  algorithms 
involve  more  complicated  numerical  operations  than  do  the  MS-E-C 
algorithms,  which  may  amplify  the  propagation  of  noise  and  errors 
into  the  reconstructed  image.  Therefore,  the  use  of  a  MS-E-C 
algorithm  instead  of  a  MS-FBPP  algorithm  (using  the  same  weight¬ 
ing  function)  will  generally  result  in  a  reduction  of  the  recon.structed 
image  variance  and/or  a  reduction  of  the  image  artifacts. 

Recently,  non-linear  reconstruction  algorithms  that  incorporate 
higher-order  scattering  approximations  have  been  proposed  for  full- 
scan  DT  (Lu  and  Zhang,  1996;  Tsihrintzis  and  Devaney,  2000b; 
Tsihrintzis  and  Devaney,  2000a).  The  generalization  of  these  works 
to  case  of  minimal-scan  DT  is  an  important  task  that  is  currently 
under  way. 
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/Ifc/rflcr— Diffraction  tomography  (DT)  is  an  inversion  technique  that 
reconstructs  the  refractive  index  distribution  of  a  scattering  object.  We  pre¬ 
viously  demonstrated  that  by  exploiting  the  redundant  information  in  the 
DT  data,  the  scattering  object  could  be  exactly  reconstructed  using  mea¬ 
surements  taken  over  the  angular  range  [0,  where  tt  <  <^min  < 

3tr /2.  In  this  paper,  we  reveal  a  relationship  between  the  maximum  scan¬ 
ning  angle  and  image  resolution  when  a  filtered  backpropagation  (FBPP) 
reconstruction  algorithm  is  employed  for  image  reconstruction.  Based  on 
this  observation,  we  develop  short-scan  FBPP  algorithms  that  reconstruct 
a  low-pass  filtered  scattering  object  from  measurements  acquired  over  the 
angular  range  (0,  where 

Index  Terms— Diffraction  tomography,  limited-view  tomography,  wave- 
field  inversion  techniques. 


I.  Introduction 

In  diffraction  tomography  (DT),  a  semi-transparent  scattering  ob¬ 
ject  is  inteiTogated  using  a  diffracting  optical  or  acoustical  wavefield 
and  the  scattered  wavefield  around  the  object  is  measured  and  used  to 
reconstruct  the  (low-pass  filtered)  refractive  index  distribution  of  the 
scattering  object.  The  principles  of  DT  have  been  extensively  utilized 
for  developing  optical  and  acoustic  tomographic  imaging  systems.  Re¬ 
cently,  interest  in  DT  within  the  optical  imaging  community  has  in¬ 
creased  because  of  its  potential  application  to  the  diffuse-photon  den¬ 
sity  wave  tomography  [l]-[3]. 

It  was  shown  previously  [4],  [5]  that,  in  two-dimensional  (2-D)  DT 
employing  plane- wave  or  cylindrical -wave  sources  and  the  classical 
scanning  geometry,  one  can  reconstruct  the  scattering  object  from  a 
minimal-scan  data  set  comprised  of  measurements  acquired  over  the 
angular  range  [0,  Omin],  where  tt  <  (puun  <  Zn/2  is  specified  by 
the  measurement  geometry.  In  this  paper,  we  reveal  a  relationship  be¬ 
tween  image  resolution  and  maximum  scan  angle,  based  upon  which 
short-scan  algorithms  can  be  designed  for  reconstructing  a  low-pass 
filtered  scattering  object  from  measurements  acquired  over  the  angular 
range  [0,  where  When  the  scattering  object  is  suf¬ 

ficiently  bandlimited,  it  can  be  exactly  reconstructed  from  the  lim¬ 
ited-view  measurements  in  [0,  We  present  numerical  examples 
that  confirm  our  theoretical  assertions. 

II.  Background 

Consider  the  classical  scanning  geometry  of  DT  with  a  cylin¬ 
drical  wave  source,  as  shown  in  Fig.  1.  Let  {x^y)  and  (r,6>) 
denote  the  fixed  Cartesian  and  polar  coordinate  systems  and 
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Measured  data  /' 


Fig.  1.  The  fan-beam  scanning  geometry  of  2-D  DT.  The  interrogating 
cylindrical  wave  propagates  along  the  y  axis  and  the  scattered  wave  field  is 
measured  along  the  line  y  =  1.  5  (or  D)  denotes  the  distance  between  the 
source  (or  the  detector)  and  the  center  of  rotation. 


(f ,  r/)  the  rotated  coordinate  system.  These  systems  are  related  by 
X  =  rcosd,  y  —  rsinO,  ^  —  xcosd  -f-  ^siiid  =  rcos(<&  -  8),  and 
r)  =  -xsincj)  -h  ycostp  =  -rsin(0  -  0).  The  scattering  object,  which 
is  embedded  in  a  lossless  and  homogeneous  background  medium, 
is  illuminated  by  a  monochromatic  cylindrical-wave  Ui(^,^)  with 
complex  amplitude  Uo  and  wavenumber  k  =  2nvo,  generated  by 
a  line  source  located  at  the  position  y  =  —  S  on  the  rj  axis.  From 
measurements  of  the  scattered  wavefield  on  the  ^  axis  at  different 
view  angles  one  seeks  to  reconstruct  the  scattering  object  function 
a(f),  which  is  related  to  the  refractive  index  distribution  n(f)  within 
the  scattering  object  by  a{f)  =  -  1. 

Let  u.(C  0)  and  <t>)  denote  the  total  and 

scattered  wavefields  measured  along  the  line  r]  =  D  oriented  at  angle 
6,  as  shown  in  Fig.  1 .  For  the  sake  of  convenience,  we  introduce  a 
modified  data  function  M [vm,<l>)  that  can  be  obtained  readily  from 
the  scattered  wavefield  and  is  defined  as 

M(u„.d>)  =  [-j2rr(i/'  -  po)D] 

where  ,\'  =  s/S/{S+  D),  v'  =  s/vl-  vljx^,  and  = 

l/(27r)  The  special  case  of  plane-wave  illumi¬ 

nation  (5  ^  oc)  corresponds  to  x  =  1.  Under  the  Bom  and  paraxial 
approximations,  Devaney  derived  the  fan-beam  Fourier  diffraction  pro¬ 
jection  (FDP)  theorem  [6],  which  relates  u(tO  to  the  modified  data 


function  by 

M  =  j 

r  a{f) 

J — oc 

•  exp  |-j27r 

v'  -  r-o)  1  df, 

if  km|  <X^O 
=0 

if  |l/m|  >  Xl'O. 

(2) 
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The  FDP  theorem  can  also  be  derived  by  employing  the  Rytov 
approximation.  In  this  case,  (2)  remains  unchanged  and  only  (1)  needs 
to  be  appropriately  redefined  [6].  When  d  is  varied  from  0  to  27r,  the 
FDP  theorem  specifies  a  circular  disk  of  (double)  coverage  centered 
at  the  origin  with  radius  I'o  \/l  +  1/x^  in  the  2*D  Fourier  space  of 
a(f).  Conventional  full-scan  reconstruction  algorithms,  such  as  the 
well  known  filtered  backpropagation  (FBPP)  algorithm  [7],  utilize  this 
Fourier  space  coverage  for  reconstmcting  a(f).  We  will  refer  to  such 
a  (low-pass  filtered)  reconstructed  a(f)  as  the  “exact”  image. 

According  to  the  fan-beam  FDP  theorem  in  (2),  the  modified  data 
function  satisfies  the  consistency  condition  [4] 

M  (i/r„ ,  0)  =  M  (-I^m >  0  +  TT  -  2a)  (3) 

where  sin  a  =  sgn(i^m)  —  i^o)^/(i^m/x^)  +  ”  ^o)  )]  . 

Using  (3),  one  can  show  [4],  [5]  that  the  minimal-scan  data  acquired 
in  the  angular  range  [0,  ^min]  specifies  a  circular  disk  (with  radius 
1^0  \/l  +  l/x^)  of  coverage  in  the  Fourier  space  of  a(f),  where 

(D,„in  =  TT  +  26  and  sin6  =  .  (4) 

III.  A  Limited-View  Reconstruction  Problem  for  2-D  DT 

We  focus  now  on  a  limited-view  problem,  in  which  data  are  acquired 
only  over  the  angular  range  [0,  where  tt  <  $‘^  <  0m:n .  In  this 
situation,  it  is  well  known  that  the  exact  image  cannot,  in  general,  be  re- 
constmcted  [8].  However,  we  demonstrate  that  algorithms  can  be  devel¬ 
oped  for  reconstructing  a  low-pass  filtered  approximation  of  the  exact 
image.  Consider  a  scattering  object  (f)  whose  2-D  Fourier  transform 
A'iu)  isbandlimitedtoadisk  ofradiusi?c  centered  at  the  origin,  where 


and  0  <  Uc  <  xi^q.  Then,  according  to  the  fan-beam  FDP  theorem 
in  (4),  the  modified  data  function  is  nonzero  only  for 

li'ml  <  I'c  The  data  space  Wc  =  [Wm\  <  <  0  <  2;r],  in 

which  the  modified  data  function  M{i'm,<i>)  is  defined,  can  be 
divided  into  the  four  subspaces  A,  B,  C,  and  V,  as  shown  in 
Fig.  2,  where  A  ==  ?0  <  ^  <  2Q{um)  +  2a(i/r:)], 

B  =  [|i^m|  <  2a{um)  +  2a(i/c)  <  (^  <  tt  -f  2a(i/„i)], 

C  =  [|i^m|  <  -h  2a{UrTz)  <  <l>  <  s^d 

T>  =  [|i/m|  <  ^  <  27r].  The  value  of  is  determined  by 

=  TT  -h  2a  (uc) .  (b) 

Using  (3),  it  can  be  verified  that  information  of  ^)  in  subspace 

A  is  redundant  to  that  of  Mii/m.d)  in  subspace  C.  Similarly,  infor¬ 
mation  of  M {Um ,  6)  in  subspace  B  is  redundant  to  that  of  M 
in  subspace  V.  Therefore,  in  principle,  the  modified  data  function 
Af(i^,„,c^)  is  completely  specified  by  its  values  in  the  subspaces  A 
and  B.  However,  because  the  boundary  between  the  subspaces  B  and 
C  is  a  nonlinear  function  of  Vm  and  6  and  because  each  horizontal 
line  in  Wc  corresponds  to  a  measurement  acquired  at  a  particular 
angle  6,  the  information  in  subspaces  B  and  C  cannot  in  practice  be 
determined  independently  of  each  other.  Consequently,  in  order  to 
determine  Af(i/m><^)  in  subspaces  A  and  B,  it  is  necessary  to  scan 
the  union  AUBUC  =  l\um\  <  t'c.O  <  4>  <  This  observation 
can  also  be  understood  by  examining  the  2-D  Fourier  space  coverage 
of  a(r)  that  is  obtained  by  varying  the  scanning  angle  from  0  to 
As  shown  in  Fig.  3,  although  the  disk  of  Fourier  space  coverage  with 


Fig.  2.  The  fiill-scan  data  space  =  AUBUCUV  contains  data  in  the 
angular  range  [0,  27r].  The  subspace  AUBUC‘m[0,  $"]  contains  all  of  die 
information  ncessary  for  exact  reconstruction  of  the  scattering  object  function 
whose  2-D  Fourier  transform  is  bandlimited  to  a  disk  of  radius  Rdt'c). 


Fig.  3 .  The  2-D  Fourier  space  coverage  of  the  scattering  object  that  is  obtained 
by  varying  6  from  0  to  The  disk  of  Fourier  space  coverage  with  radius 
R  =  1^0  yj\  -F  1/x^  is  incomplete,  with  the  shaded  region  denoting  the  missing 
data.  However,  the  coverage  corresponding  to  the  disk  of  radius  Rc ,  which  is 
defined  by  (5),  is  completely  specified.  In  generating  the  figure,  =  230® 
and  X  =  1  were  utilized. 

radius  i?  =  is  incomplete,  the  coverage  corresponding 

to  the  disk  of  radius  Rc(Pc)  is  completely  specified.  Therefore,  in 
order  to  exactly  reconstruct  whose  2-D  Fourier  transform 

is  bandlimited  to  a  disk  of  radius  Rc{vc),  only  measurements 
corresponding  to  view  angles  in  [0,  are  required.  Alternatively, 
for  an  arbitrary  scattering  object  a(f)  and  specified  >  tt,  one 
can  readily  reconstruct  which  is  a  low-pass  filtered  version  of 
a(f)  whose  2-D  Fourier  transform  is  bandlimited  to  the  disk  of  radius 
Rc{i^c),  where  the  value  of  the  data  cutoff  frequency  0  <  i^c  <  Xi'o 
is  determined  by  (6). 

A  plot  of  versus  vd  vq  for  plane-  and  cylindrical -wave  illumina¬ 
tion  is  shown  in  Fig.  4.  As  expected,  the  maximum  scanning  angle 
is  a  monotonically  increasing  function  of  the  data  cutoff  frequency  Vc. 
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The  nonlinear  shape  of  the  curves  indicates  the  scanning  angle  can  be 
reduced  from  Omin  (e.g.,  by  30°)  with  little  loss  of  resolution  in  the  re¬ 
constructed  image.  Also,  the  fact  that  the  plane-wave  (x  =  1)  curve  is 
everywhere  higher  than  the  cylindrical-wave  (x  =  0.5)  curve  reflects 
the  fact  that  the  angular  scanning  requirements  of  plane-beam  DT  are 
more  restrictive  than  for  fan -beam  DT  [5]. 

IV.  Short-Scan  Reconstruction  Algorithms  for  Limited-View 

DT 

Although  the  data  AuBUCin  Fig.  2  contains  all  of  the  informa¬ 
tion  necessary  for  reconstruction  oftt°(r),  subspaces  A  and  C  contain 
redundant  information  that  needs  to  be  properly  normalized  in  the  re¬ 
construction  process.  This  can  be  achieved  by  introducing  a  weighted 
modified  data  function  as  [4],  [5] 


(Um,4>)  =  W  (l^m,  M  [Urn,  <t>) 

(7) 

where  satisfies 

W  {Utri .  o)  -H  u’  (  —  I^rn  .  +  TT  “  2o)  =  1 

(8a) 

everywhere  in  the  data  space  Wc 

w{Urn,<l>)  =  1 

(8b) 

in  subspace  B  and 

w  (um  ,</>)  =  () 

(8c) 

in  the  subspace  {DU[|i/„.  I  >  i/c^O  <  <?>  <  2  tt]}.  The  image  a ‘■'(f)  can 
be  reconstructed  using  a  short-scan  FBPP  (SS-FBPP)  reconstruction 
algorithm  given  by 


„(”)(r,e)=  r  r 

Jc^)  Jurr,^-u^  ^ 

,|M'  (t/,„,6) 

X  exp  |j27r  sgn(r^m  ) 

i/^  +(i/'-  j/o)^r(:os(d- 

o-») 

V 

X  dvnxdo 

(9) 

which  reduces  to  the  fiill-scan  fan-beam  FBPP  algorithm  [5]  when 
=  27r  and  =  1/2.  Note  that  different  choice  for 

w{v,n^Q)  that  satisfy  (8),  in  effect,  specify  different  SS-FBPP 
algorithms. 


V  Numerical  Results 

To  validate  the  theoretical  results  above,  we  considered  a  numerical 
phantom  containing  two  elliptical  disks  whose  2-D  Fourier  transform 
was  approximately  bandlimited  to  a  disk  of  radius  Rc{^c  =  0.45i^o) 
[see  (5)].  Data  sets  of  simulated  scattered  fields  were  generated  using 
the  plane-wave  FDP  theorem  (i.e.,  x  =  1)  and  using  various  values  for 
We  reconstructed  images,  which  are  shown  in  Fig.  5,  from  these 
data  sets  using  the  conventional  FBPP  and  SS-FBPP  algorithms.  The 
SS-FBPP  algorithm  was  specified  by  a  weighting  function  w{um,4>) 
that  took  on  the  values  1/2, 1, 1/2,  and  0,  in  the  data  subspaces  A  H,  C, 
and  V,  respectively.  Fig.  5(a)  shows  images  reconstructed  by  use  of 
the  FBPP  algorithm  (left)  and  SS-FBPP  algorithm  (right),  using  data 
sets  corresponding  to  =  2;r  and  =  0riiin(x  =  1)  =  37r/2, 
respectively.  It  is  observed  that  both  images  appear  virtually  identical, 
reflecting  the  fact  that  both  of  these  data  sets  contain  the  complete  in¬ 
formation  about  the  scattering  object. 

Fig.  5(b)  shows  images  reconstructed  by  use  of  the  FBPP  algorithm 
(left)  and  SS-FBPP  algorithm  (right),  using  a  data  set  corresponding 
to  =  TT  -f  2<^((h45^(«  207°).  Clearly,  the  image  reconstructed 

using  the  FBPP  algorithm  is  distorted  and  contains  artifacts.  However, 
the  image  reconstructed  by  use  of  the  SS-FBPP  algorithm  appears  cor¬ 
rect  and  virtually  identical  to  the  images  shown  in  Fig.  5(a).  This  con¬ 
firms  our  assertion  that  the  SS-FBPP  algorithms,  which  utilize  in  the 
angular  range  [0,  can  exactly  reconstruct  a  scattering  object  a‘'(f) 

whose  2-D  Fourier  transform  A‘'(i?)  is  bandlimited  to  a  disk  of  radius 
Rr{j^c),  where  Uc  and  are  related  by  (6). 

Fig.  5(c)  shows  images  reconstructed  by  use  of  the  FBPP  algorithm 
(left)  and  SS-FBPP  algorithm  (right),  using  a  data  set  corresponding  to 
=  TT  -f  2a(0.25^'o)(~  195°).  Note  that  because  the  2-D  Fourier 
transform  of  a(r)  has  support  on  the  disk  or  radius  Rc{i^c  =  0.45j/o), 
the  measurements  in  the  angular  range  [0,  =  195°]  do  not  com¬ 

pletely  specify  the  scattering  object  (i.e.,  the  disk  of  coverage  in  2-D 
Fourier  space  with  radius  Rdi^r  =  0.45i/o)  will  not  be  completely 
filled  in.)  As  expected,  the  image  reconstructed  using  the  FBPP  al¬ 
gorithm  is  blurred,  distorted  and  contains  artifacts.  The  image  recon¬ 
structed  using  the  SS-FBPP  algorithm  also  appears  blurred,  but  does 
not  contain  any  noticeable  distortions  or  artifacts.  This  confirms  our 
assertion  that,  when  the  2-D  Fourier  transform  of  a  scattering  object 
a{f)  is  not  bandlimited  to  a  disk  of  radius  i?c(j^c),  the  SS-FBPP  algo¬ 
rithms  that  utilize  the  measurements  corresponding  to  view  angles  in 
[0,  (where  Uc  and  are  related  by  (6))  can  reconstruct  a  low-pass 
filtered  version  of  a{f)  whose  2-D  Fourier  transform  is  bandlimited  to 
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Fig.  5.  Images  reconstructed  using  the  FBPP  and  SS-FBPP  algorithms  for 
various  simulated  data  sets.  See  the  text  for  a  detailed  description. 


the  disk  of  radius  In  this  particular  example,  the  2-D  Fourier 

transform  of  the  image  reconstructed  by  use  of  the  SS-FBPP  algorithm 
is  bandlimited  to  a  disk  of  radius  !?<:(. 25/ .45  r^c). 

VI.  Summary 

We  demonstrated  previously  [4],  [5]  that  in  2-D  DT  employing 
plane-wave  or  cylindrical-wave  sources,  one  can  exactly  reconstruct 
the  scattering  object  from  a  minimal-scan  data  set  acquired  using  view 
angles  only  in  [0,  ^min],  where  n  <  <l>„un  <  37r/2  is  a  specified 
function  of  the  measurement  geometry.  In  this  study,  we  have  demon¬ 
strated  that  when  measurements  are  available  only  for  view  angles  in 
[0,  where  tt  <  <  ^min,  a  simple  relationship  exists  between 

the  maximum  scanning  angle  and  the  image  resolution  when  a 
FBPP  algorithm  is  employed  to  reconstruct  the  image.  By  properly 
weighting  the  measurement  data,  a  low-pass  filtered  approximation 
of  the  scattering  object  that  is  free  of  conspicuous  artifacts  can  be 
obtained  from  the  measurements  corresponding  to  view  angles  in 
[0,  When  the  scattering  object  is  sufficiently  bandlimited,  it 
can  be  exactly  reconstructed.  This  observation  is  practically  useful, 
because  it  provides  a  convenient  mechanism  for  regularizing  the 
severely  ill-posed  limited-view  DT  reconstruction  problem;  when  the 
maximum  scanning  angle  is  greater  than  tt,  a  stable  reconstruction 
can  always  be  performed  by  sacrificing  spatial  resolution  in  the  recon¬ 
structed  image.  It  can  be  demonstrated  that  the  statistical  properties 
of  the  SS-FBPP  algorithms  are  qualitatively  similar  to  those  of  the 
minimal-scan  FBPP  reconstruction  algorithms  investigated  previously 
[5].  In  the  limited-view  radon  transform  inversion  problem  [9],  an 
analogous  regularization  mechanism  does  not  exist  and  some  sort  of  a 
priori  information  regarding  the  object  function  is  generally  required 
to  effectively  regularize  the  problem. 

Because  we  have  assumed  a  2-D  imaging  geometry  in  this  study,  the 
developed  SS-FBPP  reconstruction  algorithms  may  be  useful  for  appli¬ 
cations  in  which  out-of-plane  scattering  is  not  significant.  In  diffuse- 
photon  density  wave  tomography,  the  wavenumber  is  complex-valued 
and  the  FDP  theorem  describes  a  mapping  between  the  data  function 
and  a  set  of  complex-valued  frequencies  of  the  scattering  object  func¬ 
tion’s  Fourier  transform.  The  extension  of  the  concepts  and  techniques 
introduced  in  this  correspondence  to  the  case  where  the  wavenumber 
is  complex-valued  and  to  the  three-dimensional  reconstruction  problem 
represent  important  topics  for  future  research. 
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